Phonon-assisted luminescence by finite atomic displacements-Yambopy

[Sam-2025]

Dear all,

I am also checking the following tutorial :
Phonon-assisted luminescence by finite atomic displacements
for cubic BTO, it is indirect and the q = 1/2 , 1/2 , 1/2 ,
first I run scf, nscf, phonon and got Q.mode as attached for q= 1/2 , 1/2 , 1/2

that is why I created the diagonal supercell (using the attached python script )

Q.modes.txt

, but then when I want to create the displacement structure, it is giving the following error , even though Q.mode is 225 ,

Could you please guide me , why I am getting this error ?

(yamboenv) [sghanbar@vacc-login1 PL]$ python generate_supercells.py -S
Applying displacements according to phonon modes...
Traceback (most recent call last):
File "/gpfs1/home/s/g/sghanbar/Yambo/PL/generate_supercells.py", line 142, in <module>
if args.displaced_nondiagonal: generate_displaced_supercells(uc_filnm,Q,modes_file,kpoints)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/gpfs1/home/s/g/sghanbar/Yambo/PL/generate_supercells.py", line 82, in generate_displaced_supercells
sc.displace(eivs,nd_atom_positions,Temp=0.1) # Generate list of displaced supercells as PwIn objects called 'modes_qe'
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/users/s/g/sghanbar/yamboenv/lib/python3.11/site-packages/qepy/supercell.py", line 108, in displace
self.initialize_phonons(modes_file,self.atypes,Temp)
File "/users/s/g/sghanbar/yamboenv/lib/python3.11/site-packages/qepy/supercell.py", line 176, in initialize_phonons
self.eigs = read_eig(modes_file,self.basis)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/users/s/g/sghanbar/yamboenv/lib/python3.11/site-packages/qepy/supercell.py", line 49, in read_eig
eig = np.reshape(eig, (modes,basis,3))
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/users/s/g/sghanbar/yamboenv/lib/python3.11/site-packages/numpy/_core/fromnumeric.py", line 324, in reshape
return _wrapfunc(a, 'reshape', shape, order=order)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/users/s/g/sghanbar/yamboenv/lib/python3.11/site-packages/numpy/_core/fromnumeric.py", line 57, in _wrapfunc
return bound(*args, **kwds)
^^^^^^^^^^^^^^^^^^^^
ValueError: cannot reshape array of size 222 into shape (15,5,3)

[claudio]

Dear Samira

may you put in the forum also your DFT input file, so I can reproduce the error
best
Claudio

[Sam-2025]

Dear professor,

For sure. Thank you so much.
In cubic BaTiO₃, the band gap is indirect:

The valence band maximum (VBM) is typically at the R point (½, ½, ½)

The conduction band minimum (CBM) is at the Γ point (0, 0, 0)

For BaTiO₃ (BTO) using q = (½, ½, ½) and a 2×2×2 diagonal supercell( To model phonons at q , we need a supercell that folds this q-point to Γ, In reciprocal space, (½, ½, ½) folds to Γ if the real-space lattice expands 2× in each direction, that is why we need a diagonal supercell, I thought), then I modified the python script on the tutorial to get the displaced structures for diagonal supercell, and attached for your consideration. I used "python generate_displaced_supercell.py -D" and create the supercell with 80 atoms (80 atoms is a lot and computationally expensive, but I don't know what to do ) and then "python generate_displaced_supercell.py --displaced_diag" to get the displacement structures, It must gives us 15 displaced structure, because for q point , we have 15 freq, (3N), but I don't know why I am getting 12 structure ! (to solve the error :ValueError: cannot reshape array of size 222 into shape (15,5,3), I manually add a line in matdyne.modes file) but then got 12 displaced structure not 15! which is not correct again.

my SCF file :
&control
calculation='scf',
prefix='bto',
restart_mode='from_scratch'
pseudo_dir = './Pseudos'
wf_collect=.true.
nstep = 200
verbosity = 'high'
/
&system
ibrav = 1,
celldm(1)=7.584
nat= 5,
ntyp= 3,
force_symmorphic=.true.
ecutwfc = 80,
occupations = 'smearing'
degauss=0.01
/
&electrons
electron_maxstep = 100
conv_thr = 1.0e-8
/
ATOMIC_SPECIES
Ba 137.327 Ba_std.upf
Ti 47.87 Ti_std.upf
O 15.99 O_std.upf
ATOMIC_POSITIONS {crystal}
Ba 0.00 0.00 0.00
Ti 0.50 0.50 0.50
O 0.00 0.50 0.50
O 0.50 0.50 0.00
O 0.50 0.00 0.50
K_POINTS {automatic}
6 6 6 0 0 0

nscf file :
&control
calculation='nscf',
prefix='bto',
pseudo_dir = './Pseudos'
nstep = 200
verbosity = 'high'
wf_collect=.true.
/
&system
ibrav = 1,
celldm(1)=7.584
nat= 5,
ntyp= 3,
force_symmorphic=.true.
ecutwfc = 80,nbnd = 80
occupations = 'smearing'
degauss=0.01
/
&electrons
diago_thr_init=5.0e-6,
diago_full_acc=.true.
electron_maxstep = 100,
diagonalization='david'
mixing_mode = 'plain'
mixing_beta = 0.6
conv_thr = 1.0e-8
/
ATOMIC_SPECIES
Ba 137.327 Ba_std.upf
Ti 47.87 Ti_std.upf
O 15.99 O_std.upf
ATOMIC_POSITIONS {crystal}
Ba 0.00 0.00 0.00
Ti 0.50 0.50 0.50
O 0.00 0.50 0.50
O 0.50 0.50 0.00
O 0.50 0.00 0.50
K_POINTS {automatic}
6 6 6 0 0 0

phonon.in
bn_phonons
&inputph
tr2_ph=1.0d-13,
verbosity='high'
prefix='bto',
outdir='./',
fildyn='bto.dynG'
ldisp=.true.
nq1=3, nq2=3, nq3=3
/

q2r.in

&input
fildyn='bto.dynG'
zasr='crystal',
flfrc='bto.fc'
/


matdyne.modes :
&input
asr='crystal',
flfrc='bto.fc',
flfrq='bto.freq',
deltaE=1.d0,
q_in_cryst_coord=.true.
fleig='bto.eig'
/
1
0.5 0.5 0.5


I also attached the matdyne.modes files for your consideration. Professor, I sincerely appreciate it.
Professor, I add the following in the python script to generate the displaced structure for diagonal supercell :

def generate_displaced_diagonal_supercell(uc,R,eivs):
"""
New: displaced structures from diagonal supercell
"""
qe = PwIn.from_file(uc)
sc = Supercell(qe)
atom_positions = sc.d_sup(R) # Make diagonal supercell

sc.displace(eivs, atom_positions, Temp=0.1)
N_modes = len(sc.modes_qe)

suffix = ''.join(sc.qe_d.control['calculation'].split('\''))
prefix = ''.join(sc.qe_d.control['prefix'].split('\''))

for mode in range(N_modes):
sc.modes_qe[mode].write('sc_diag_displaced_%s_mode%d.%s'%(prefix, mode+1, suffix))

print('Displaced diagonal supercells written to file.')

and then :
parser.add_argument('--displaced_diag', action="store_true", help='Build displaced diagonal supercell')

finally :
if args.displaced_diag: generate_displaced_diagonal_supercell(uc_filnm, R, modes_file)

[Sam-2025]

Dear Professor,

So sorry to bother you a lot. For this tutorial:https://wiki.yambo-code.eu/wiki/index.p ... ite_note-7, if I got the displaced structures, I do not need to run again scf calculation? only need to do an nscf ? and then bse.in ? Thank you so much for your guidance.


Professor, I also used the attached python script and also in my yamboenv, I changed the restriction to create the displaced structure based on the freq, and it gave me 15 displaced structures. (I patched qepy’s supercell.py to remove its built-in mode filtering (specifically deleting the if GAMMA: self.disps = self.disps[3:] line), ensuring none of the 15 R-point phonon modes in cubic BaTiO₃ are skipped. As a result, qepy now generates all 15 displaced supercells for use in Yambo’s finite-displacement luminescence workflow.) I am not sure how much this way is correct. Could please also look at this one? Thank you so much.

# -*- coding: cp1252 -*-
import numpy as np
np.float_ = np.float64

from qepy import *
#from yambopy.lattice import *
import sys
import argparse
"""
In this example we construct supercells (sc) of bulk hBN.
Starting from a unit cell (uc) input that is read from file, we build:

- diagonal supercell
- non-diagonal supercell folding a specific q-point Q
- displaced non-diagonal supercells along phonon eigenmodes at Q (read from file)
"""

def generate_diagonal_supercell(uc,R):
"""
First case: diagonal supercell of size R
"""
qe = PwIn.from_file(uc) #read the uc input

sc = Supercell(qe) # initialize class
sc.d_sup(R) # Generate supercell as PwIn object called 'qe_d'

#name of output file
suffix = ''.join(sc.qe_d.control['calculation'].split('\''))
prefix = ''.join(sc.qe_d.control['prefix'].split('\''))

#write supercell to file
sc.qe_d.write('sc_diagonal_%s.%s'%(prefix,suffix))

print('Diagonal supercell written to file.')

def generate_nondiagonal_supercell(uc,Q,kpoints=None):
"""
Second case: nondiagonal supercell folding point Q
"""
qe = PwIn.from_file(uc) #read the uc input
if kpoints is not None: qe.kpoints = kpoints #Optional: manually change the original kpt mesh to ensure consistency

sc = Supercell(qe) # initialize class
sc.nd_sup(Q) # Generate supercell as PwIn object called 'qe_nd'

#name of output file
suffix = ''.join(sc.qe_nd.control['calculation'].split('\''))
prefix = ''.join(sc.qe_nd.control['prefix'].split('\''))

#write supercell to file
sc.qe_nd.write('sc_nondiagonal_%s.%s'%(prefix,suffix))

print('Nondiagonal supercell written to file.')


def generate_displaced_Rpoint_supercell(uc, eivs):
"""
Fifth case: displaced 2×2×2 diagonal supercell for q = (1/2,1/2,1/2) (R-point)
"""
qe = PwIn.from_file(uc) # read the unit cell pw input

sc = Supercell(qe) # initialize Supercell
atom_positions = sc.d_sup([2,2,2]) # build 2×2×2 supercell (diagonal)

# Displace atoms using phonon eigenmodes
sc.displace(eivs, atom_positions, Temp=0.1)
N_modes = len(sc.modes_qe)

suffix = ''.join(sc.qe_d.control['calculation'].split('\''))
prefix = ''.join(sc.qe_d.control['prefix'].split('\''))

# write displaced supercells to file
for mode in range(N_modes):
sc.modes_qe[mode].write('sc_R_displaced_%s_mode%d.%s' % (prefix, mode+1, suffix))

print('Displaced R-point supercells written to file.')



def generate_displaced_supercells(uc,Q,eivs,kpoints=None):
"""
Third case: displaced supercells along phonon modes at Q
"""
qe = PwIn.from_file(uc) #read the uc input
if kpoints is not None: qe.kpoints = kpoints #Optional: manually change the original kpt mesh to ensure consistency

sc = Supercell(qe) # initialize class
nd_atom_positions = sc.nd_sup(Q) # Generate supercell as PwIn object called 'qe_nd' getting new atomic positions

# Displace atoms.
# Intensity is Temp (in bohr)
# Sign and direction (standing wave at Q) given by modes_file
sc.displace(eivs,nd_atom_positions,Temp=0.1) # Generate list of displaced supercells as PwIn objects called 'modes_qe'
N_modes = len(sc.modes_qe)

#name of output file
suffix = ''.join(sc.qe_nd.control['calculation'].split('\''))
prefix = ''.join(sc.qe_nd.control['prefix'].split('\''))

#write supercells to file for each mode
for mode in range(N_modes): sc.modes_qe[mode].write('sc_displaced_%s_mode%d.%s'%(prefix,mode+1,suffix))

print('Displaced supercells written to file.')

def generate_displaced_unitcell(uc,eivs):
"""
Fourth case: displaced cell at Q=0
"""
qe = PwIn.from_file(uc) #read the uc input

sc = Supercell(qe) # initialize class
atom_positions = sc.d_sup([1,1,1]) # Generate "supercell" with size 1 as PwIn object called 'qe_d' getting atomic positions

# Displace atoms.
# Intensity is Temp (in bohr)
# Sign and direction (standing wave at Q) given by modes_file
sc.displace(eivs,atom_positions,Temp=0.1) # Generate list of displaced supercells as PwIn objects called 'modes_qe'
N_modes = len(sc.modes_qe)

#name of output file
suffix = ''.join(sc.qe_d.control['calculation'].split('\''))
prefix = ''.join(sc.qe_d.control['prefix'].split('\''))

#write supercells to file for each mode
for mode in range(N_modes): sc.modes_qe[mode].write('uc_GAMMA_displaced_%s_mode%d.%s'%(prefix,mode+4,suffix))

if __name__ == "__main__":
parser = argparse.ArgumentParser(description='Supercell generation')
parser.add_argument('-D', '--diagonal', action="store_true", help='Build diagonal supercell')
parser.add_argument('-N', '--nondiagonal', action="store_true", help='Build nondiagonal supercell')
parser.add_argument('-R', '--displaced_Rpoint', action="store_true", help='Build displaced supercell for q=(1/2,1/2,1/2) (R-point)')
parser.add_argument('-S', '--displaced_nondiagonal', action="store_true", help='Build displaced nondiagonal supercells')
parser.add_argument('-G', '--displaced_gamma', action="store_true", help='Build displaced unit cell')
args = parser.parse_args()

if len(sys.argv)==1:
parser.print_help()
sys.exit(1)

uc_filnm = 'bto_nscf.in' #pw input
R = [2,2,2] #diagonal supercell dimensions
Q = [[1,1,1],[2,2,2]] #qpoint to fold in nondiagonal supercell: Q=(1/3,-1/6,0)
kpoints = [6,6,6] #NB: these fractional crystal coordinates must exactly divide the kpoint mesh!
modes_file = 'R.modes' #file with phonon eigenmodes at Q. NB2: we use those already scaled with the atomic masses

if args.diagonal: generate_diagonal_supercell(uc_filnm,R)

if args.nondiagonal: generate_nondiagonal_supercell(uc_filnm,Q,kpoints)

if args.displaced_Rpoint:
generate_displaced_Rpoint_supercell(uc_filnm, modes_file)

#NB3: If computing the eigenmodes at Q with qe, remember to use Q in *Cartesian* coordinates
# in both qe input and output!!
# For example, Q_crystal=(1/3,-1/6,0) ====> Q_cartesian = (1/3,0,0) in the BN hexagonal lattice

if args.displaced_nondiagonal: generate_displaced_supercells(uc_filnm,Q,modes_file,kpoints)

#NB4: For the following we are using the same modes_file as in the case above for simplicity, but it is *wrong*:
# Phonon modes at Gamma must be used in this case.
if args.displaced_gamma: generate_displaced_unitcell(uc_filnm,modes_file)

[claudio]

Dear Samira

yes, you should generate 15 displaced cell.
Now you can run SCF + NSCF in all supercells.
In the undisplaced one you run the calculation of dielectric constant + BSE.

Then you can copy the dielectric constant, files SAVE/ndb.em1s\* in the other suprcells and run the BSE

best
Claudio