Question about ndb.dipole data structure with open shell system

[Shimin]

Dear Yambo developers,

I'm trying to get the transition dipole moment <i|[r,H]|j> from the the ndb.dipoles

I have a question regarding the each dimension of the data. Here i show the datashape of the dipole.

Code: Select all

database = Dataset("./ndb.dipoles","r")
print(database.variables["DIP_P"])

Code: Select all

>> float32 DIP_P(D_0000000002, D_0000000001, D_0000001242, D_0000001760, D_0000000003, D_0000000002)
unlimited dimensions:
current shape = (2, 1, 1242, 1760, 3, 2)
filling on, default _FillValue of 9.969209968386869e+36 used

If i understand correctly (2, 1, 1242, 1760, 3, 2) are corresponded to (spin up/dn , k points, v, c, direction x y z, Real and Imaginary component)

My question is when computing the spin-polarized calculation we have different numbers of valence band v and conduction band c at two spin channel.
for example my system is a open shell system that have 1242 electrons at spin up channel and 1240 electrons at spin down channel.

It is unclear to me how to get, for example <v|p|c> = <band 1240|p|band 1241> at spin down channel with current data shape. because both v and c of this matrix element are defined as v in the data.

[Shimin]

Hi,
Does anyone have an idea on this problem?

[Daniele Varsano]

Dear Shimin,

you can have a look at the Yambo-py tutorial where it is explained how to extract dipoles element from the ndb.dipoles database.

https://www.yambo-code.eu/wiki/index.ph ... absorption

Best,
Daniele

[Shimin]

Dear Daniele,

Thank you for the reply!

I was originally trying to use yambopy to extract the dipole,
but i come up with error a the first step:

Code: Select all

ylat = YamboLatticeDB.from_db_file(filename='./ns.db1')

>>

Code: Select all

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[9], line 1
----> 1 ylat = YamboLatticeDB.from_db_file(filename='./ns.db1')

File ~/Python_Env/west/lib/python3.11/site-packages/yambopy/dbs/latticedb.py:52, in YamboLatticeDB.from_db_file(cls, filename, Expand, expand_mode, atol)
     50     flatten = lambda l: [item for sublist in l for item in sublist]
     51     atomic_numbers = flatten([[tmp_an[n]]*na for n,na in enumerate(natoms_a)])
---> 52     atomic_positions = np.vstack([[tmp_apos[n,ia] for ia in range(na)] for n,na in enumerate(natoms_a) ])
     54     args = dict( atomic_numbers       = atomic_numbers,
     55                  car_atomic_positions = atomic_positions,
     56                  sym_car              = database.variables['SYMMETRY'][:],
   (...)
     59                  alat                 = database.variables['LATTICE_PARAMETER'][:].T,
     60                  time_rev             = dimensions[9] )
     62 y = cls(**args)

File ~/Python_Env/west/lib/python3.11/site-packages/numpy/core/shape_base.py:289, in vstack(tup, dtype, casting)
    287 if not isinstance(arrs, list):
    288     arrs = [arrs]
--> 289 return _nx.concatenate(arrs, 0, dtype=dtype, casting=casting)

ValueError: all the input array dimensions except for the concatenation axis must match exactly, but along dimension 1, the array at index 0 has size 0 and the array at index 1 has size 3

So i used the method readDB written in "dipolesdb.py", where it gives the dipole matrix element in following index:
[nspin, nkpoints, cartesian directions, nbands conduction, nbands valence]

where the dimension of data is :

Code: Select all

(2, 1, 3, 1760, 1242)

It still remains puzzle to me that if i have vbm=1242 at spin up channel, and vbm=1240 at spin dn channel. How can i get the transition dipole of <v|p|c> = <1240|p|1241> at spin down?

[palful]

Dear Shimin,

It may be that you are the first person using yambopy with an open-shell system!

May I ask you for some minimal test attachments, so I can have a look?

First of all, the initial error you get while opening the ns.db1 with yambopy is very strange: would you mind attaching your ns.db1 database so that I can try to check what's going on? If it is too large, maybe you could have a smaller system where you still have the same problem while loading it?

Same goes for ndb.dipoles: do you think you could send me the ns.db1 and ndb.dipoles of the smallest open-shell test system that you can do? For example, something with about 12 spin-up and 10 spin-down valence bands? (up to around 102 and 100 would also be fine, but not thousands)

Thank you!

Cheers,
Fulvio

[Shimin]

Dear Fulvio,

Thank you ! I attached my database file and scf.in in this google drive: https://drive.google.com/drive/folders/ ... sp=sharing

Best,
Shimin

[palful]

Dear Shimin,

So I checked and I think yambo is correctly computing the transition dipoles for both spin channels, even though in the database they are stored with the dimensions of the spin=0 case (i.e., 1242 valence bands instead of 1240 also for the spin=1 case). So for example the transition <1240|d|1241> should be correctly computed.

In order to show this, I have patched yambopy so that it can now read correctly the dipoles for an open-shell system (it was not supported before). In order to get it, you need to clone it from git following the "Installation from latest patch" option on the wiki.

Let me discuss more in detail below in order to fully answer your question.

- The error when loading the databases with yambopy was due to the fact that in your calculation you had set more atomic types than you actually used: in particular, you had the "Zn1" type defined but not actually used in the atomic positions. Now yambopy can handle this case and opens the database correctly. Of course, the cleanest way would be to have the correct ntyp value in your input.

- For the dipoles, I patched the open-shell case. Yambopy by defaults restructures the dipoles array into a hermitian matrix Nbands x Nbands, where the diagonal blocks (corresponding to c->c' and v->v' transitions) are zero by construction, while the off-diagonal blocks contain the v->c transition amplitudes.

In the case of your system, if I plot a slice of the dipole matrix from states 1236 to 1246 for both spins I obtain the results below (showing the norm of the dipoles with only two significant digits):

Code: Select all

# Dipoles
# SPIN 0
#      1236  1237  1238  1239  1240  1241  1242    1243  1244  1245  1246
# 1236 0.00  0.00  0.00  0.00  0.00  0.00  0.00    0.04  0.13  0.34  0.73
# 1237 0.00  0.00  0.00  0.00  0.00  0.00  0.00    0.15  0.15  2.28  0.05
# 1238 0.00  0.00  0.00  0.00  0.00  0.00  0.00    0.01  0.09  2.11  0.65
# 1239 0.00  0.00  0.00  0.00  0.00  0.00  0.00    0.09  0.01  2.11  0.65
# 1240 0.00  0.00  0.00  0.00  0.00  0.00  0.00    0.59  0.59  0.49  0.86
# 1241 0.00  0.00  0.00  0.00  0.00  0.00  0.00    0.00  0.54  0.55  0.02
# 1242 0.00  0.00  0.00  0.00  0.00  0.00  0.00    0.54  0.00  0.55  0.02

# 1243 0.04  0.15  0.01  0.09  0.59  0.00  0.54    0.00  0.00  0.00  0.00
# 1244 0.13  0.15  0.09  0.01  0.59  0.54  0.00    0.00  0.00  0.00  0.00
# 1245 0.34  2.28  2.11  2.11  0.49  0.55  0.55    0.00  0.00  0.00  0.00
# 1246 0.73  0.05  0.65  0.65  0.86  0.02  0.02    0.00  0.00  0.00  0.00

# SPIN 1
#      1236  1237  1238  1239  1240    1241  1242  1243  1244  1245  1246
# 1236 0.00  0.00  0.00  0.00  0.00    0.35  0.52  1.37  0.02  0.01  0.03
# 1237 0.00  0.00  0.00  0.00  0.00    2.28  0.00  0.13  0.09  0.09  0.01
# 1238 0.00  0.00  0.00  0.00  0.00    2.09  0.67  0.06  0.02  0.03  0.08
# 1239 0.00  0.00  0.00  0.00  0.00    2.09  0.67  0.06  0.03  0.02  0.08
# 1240 0.00  0.00  0.00  0.00  0.00    0.34  0.84  0.13  1.10  1.10  0.12

# 1241 0.35  2.28  2.09  2.09  0.34    0.00  0.00  0.00  0.00  0.00  0.00
# 1242 0.52  0.00  0.67  0.67  0.84    0.00  0.00  0.00  0.00  0.00  0.00
# 1243 1.37  0.13  0.06  0.06  0.13    0.00  0.00  0.00  0.00  0.00  0.00
# 1244 0.02  0.09  0.02  0.03  1.10    0.00  0.00  0.00  0.00  0.00  0.00
# 1245 0.01  0.09  0.03  0.02  1.10    0.00  0.00  0.00  0.00  0.00  0.00
# 1246 0.03  0.01  0.08  0.08  0.12    0.00  0.00  0.00  0.00  0.00  0.00 

As you can see, in the SPIN-0 case we see that the occupied/valence block ends correctly at 1242, while in the SPIN-1 case it ends at 1240. In particular, for the SPIN-1 case, we have nonzero values for the transition dipole |<1240|d|1241>|=0.34 (and more in general the rows 1241 and 1242 are nonzero). This for me shows that everything is working fine, but of course please check carefully that the results are what you expect.

I attach the script I used to get this file just in case you need it.

Cheers,
Fulvio