DFPT phonon calculation in monolayers

Hi everyone,

I am trying to do phonon calculation in a STO monolayer. I initially relaxed the cell by fixing the direction along the vacuum, optcell=9 and used the relaxed coordinates for the phonon calculations.

This is the input script I used for GM phonon:

  ndtset  3

#Ground state calculation
  kptopt1   1     
  tolvrs1   1.0d-10   

#Response Function calculation : d/dk
  rfelfd2   2     
   rfdir2   1 1 0   
    nqpt2   1
     qpt2   0.0 0.0 0.0  
  getwfk2   -1    
  kptopt2   2   
    iscf2  -3    
  tolwfr2   1.0d-8   
   
#Response Function calculation : electric field perturbation and phonons
  rfphon3   1           
  rfelfd3   3          
   rfdir3   1 1 0    
    nqpt3   1
     qpt3   0.0 0.0 0.0  
  getwfk3   -2      
  getddk3   -1          
  kptopt3   2        
  tolvrs3   1.0d-8

acell    1.0045392115E+01  1.0045392115E+01  4.4196407571E+01

rprim    1.0000000000E+00  0.0000000000E+00  0.0000000000E+00
         0.0000000000E+00  1.0000000000E+00  0.0000000000E+00
         0.0000000000E+00  0.0000000000E+00  1.0000000000E+00
        
natom    14 
ntypat   3   
typat    1 1 1 1 2 2 3 3 3 3 3 3 3 3   
znucl    38 22 8   
xred     6.0703843972E-21  5.0000000000E-01  4.1807979440E-01
                       5.0000000000E-01  7.4718483788E-21  4.1807979440E-01
                       5.0000000000E-01  7.4718483821E-21  5.8192020560E-01
                       6.0703844039E-21  5.0000000000E-01  5.8192020560E-01
                       5.0000000000E-01  5.0000000000E-01  5.0000000000E-01
                       6.0703843981E-21  7.4718483817E-21  5.0000000000E-01
                       8.0073670562E-01  3.0073670562E-01  5.0000000000E-01
                       6.9926329438E-01  8.0073670562E-01  5.0000000000E-01
                       3.0073670562E-01  1.9926329438E-01  5.0000000000E-01
                       1.9926329438E-01  6.9926329438E-01  5.0000000000E-01
                       5.0000000000E-01  5.0000000000E-01  4.1494586031E-01
                       6.0703843959E-21  7.4718483813E-21  4.1494586031E-01
                       6.0703843955E-21  7.4718483821E-21  5.8505413969E-01
                       5.0000000000E-01  5.0000000000E-01  5.8505413969E-01

pp_dirpath "~/pseudo/PBEsol/"
pseudos "Sr.psp8, Ti.psp8, O.psp8"

ecut     30.0 
ngkpt    8 8 1    
kptopt   1
nshiftk  1
shiftk   0.0 0.0 0.0

iprcel 45
fband 0.2

nstep    500        
diemac   6.0

however there is convergence issue with this script for the tolwfr not reaching the tolerance set even after 500 steps,
and thus giving wrong phonon frequency

scprqt: WARNING -
nstep= 500 was not enough SCF cycles to converge;
maximum residual= 8.427E-02 exceeds tolwfr= 1.000E-08

Phonon frequencies in cm-1    :
- -3.339734E-02 -2.930099E-02  0.000000E+00  0.000000E+00  0.000000E+00
-  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
-  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
-  0.000000E+00  9.162324E+01  9.162515E+01  1.047970E+02  1.048074E+02
-  1.421056E+02  1.421067E+02  1.489311E+02  1.489398E+02  2.117343E+02
-  2.117416E+02  2.316880E+02  2.317312E+02  2.492741E+02  2.897962E+02
-  2.898033E+02  4.777024E+02  4.817118E+02  4.817292E+02  5.077103E+02
-  5.077118E+02  5.174063E+02  5.174078E+02  5.751804E+02  5.751841E+02
-  5.984614E+02  8.146773E+02

chkph3 : WARNING -
Dynamical matrix incomplete, phonon frequencies may be wrong.
Likely due to a list of perturbations, as defined by rfatpol and rfdir, that does not include all displacements of all atoms and (if non-metallic material) electric field type perturbation.
Then, the dynamical matrix includes zeroes when the matrix element is not computed.
This is allowed for testing purposes. But the phonon frequencies may be wrong.
If there are symmetries, perhaps these matrix elements are zero by symmetry anyhow, and phonon frequencies might be right.
Please check the input variables rfatpol and rfdir, to determine whether abinit is doing what you intend it to do.
Note that ANADDB is able to detect whether the symmetries allow one to reconstruct the full dynamical matrix from an incomplete one. In this case, passing to ANADDB the delivered _DDB file might confirm (or not) that phonon frequencies are right.

How should such calculations be handled?

Thanks
Dominic

Hi Dominik,

we need a bit more detail. I presume it’s DS3 that stalls. Why are you running only rfdir 1 1 0? You should set 1 1 1, and get the full phonon set. It will use symmetry automatically, and z is a distinct axis wrt x and y, so you’ll definitely need it. What spgroup does it find?

You can also try to add a few more empty bands, that will help, and increase nbdbuf to have it ignore the last few bands in comparing to tolwfr. From your post we can not tell if the problem is valence bands not converging or empty (unimportant) conduction ones.

If your structure is not sufficiently well relaxed, or if it’s ferroelectric and stuck in a higher symmetry structure, the unstable phonon mode may be correct.

To make the scf converge better you can try the same preconditioning tricks as in the GS, with diemac 10, and/or reducing diemix.

Thank you for the helpful suggestions.

For the structural relaxation, I initially used optcell = 9 to relax only the in-plane lattice parameters while keeping the large vacuum layer fixed along the z-direction. I then used the relaxed coordinates obtained from this calculation for the phonon calculations. I have attached the relaxation input and output files for reference.

I also tested the calculations using rfdir 1 1 1, but I am still encountering convergence issues in DS2, where the tolwfr criterion is not reached. In this case, the phonon frequencies appear unphysical (on the order of −9000 cm⁻¹), which I suspect is due to the non-converged wavefunctions.

Regarding symmetry, the bulk structure had the I4/mcm space group, whereas the monolayer structure now shows P4/mbm, likely due to symmetry breaking along the z-direction.

Currently, I am running additional calculations by increasing the number of empty bands (testing up to 8 empty bands along with nbdbuf = 8) and using diemac=10 to check whether this improves the tolwfr convergence behavior.

Please let me know if you think there are other parameters or strategies I should consider testing.

latt_opt.abi (1.9 KB)
latt_opt.abo (372.8 KB)
phonon.abi (4.1 KB)
phonon.abo (479.0 KB)

Thanks
Dominic

ok, it was not clear this was for DS2 (tolwfr is a legal option in all cases).

This lack of convergence is quite common in DDK calculations, see the rest of the forum for other examples. The problem is that the ddk (or a perturbative E field) can in principle populate some of the conduction bands, even if you have a large gap. Imagine a large supercell with a very small E field: if you go far enough, the E.distance can be equal to the band gap. In the perturbative limit E→0, but the convergence is still complicated by admixion of conduction band character in the DDK wfk.

The solution is to add a few empty bands, and use a metallic occupation (occopt 7 for example, with tsmear 0.001). In my experience this always makes the convergence very smooth.

Note that after this, certain DFPT calculations can complain if they expect only valence bands, no unoccupied states or occopt. You can later reduce nband to nvalence, and just chop off the DDK states which are not needed (I always keep everything).

I added empty bands and metallic smearing but still now the phonon frequencies are showing -6000 cm-1. The calculations have converged for the corresponding tolerance set. Is there any mistake in the runs? I am attaching the relaxation and the phonon calculations files here.

phonon_new.abo (549.6 KB)
phonon_new.abi (4.1 KB)
latt_opt.abo (489.1 KB)
latt_opt.abi (1.9 KB)

I was able to get proper phonon frequencies with increasing the number of empty bands. The occupied states were 56 and I set nbands to 70 with nbdbuf set to 14.

I did not had to use metallic smearing in this case. The tolerance for convergence had to be lowered to 1e-9. The current GM phonons show the system is stable.

Phonon frequencies in cm-1    :
- -1.821167E-02 -1.816082E-02  0.000000E+00  7.736612E+01  7.736612E+01
-  1.035017E+02  1.035017E+02  1.053608E+02  1.175114E+02  1.175114E+02
-  1.462536E+02  1.468992E+02  1.468992E+02  1.530431E+02  2.015822E+02
-  2.015822E+02  2.085585E+02  2.085585E+02  2.267095E+02  2.267095E+02
-  2.297280E+02  2.521151E+02  2.571775E+02  2.571775E+02  2.572735E+02
-  2.853566E+02  2.853566E+02  3.465814E+02  3.524438E+02  3.883131E+02
-  4.258333E+02  4.258333E+02  4.581918E+02  4.777105E+02  4.777105E+02
-  4.779581E+02  5.199602E+02  5.677841E+02  5.677841E+02  5.679053E+02
-  5.910676E+02  8.250477E+02

Is it possible to know these calculations are correct?

I am attaching my files here for reference.

Thanks
Dominic
latt_opt.abi (2.0 KB)
latt_opt.abo (489.3 KB)

phonon.abo (407.5 KB)
phonon.abi (4.1 KB)

Hi,

I was doing finite-q phonon calculations in a STO monolayer. However I am not able to converge the phonon calculations to even 5e-10 for the tolvrs value even after 500 steps. Is there a better way rather than to increase the number of steps even further?

Am I missing some input flag that could potentially solve this convergence issue?

#Multi-dataset mode
ndtset  2

#Ground state calculation
#######################################################################
kptopt1 1               # Automatic generation of k points, taking into account the symmetry
tolvrs1 1.0d-10       # SCF stopping criterion

#Response-function calculation, phonons
#######################################################################
nqpt2   1               # If nqpt = 1, one is not doing a ground-state calculation
                        # One wavevector is to be considered
rfphon2   1             # Will consider phonon-type perturbation
qpt2   1./16. 1./16. 0.0      # This is a calculation at the Gamma point

getwfk2   -1            # Uses as input the output wf of the previous dataset

kptopt2   3             
tolvrs2 5.0d-10  # SCF stopping criterion

acell      1.0104368790E+01  1.0104368790E+01  4.0000000000E+01 Bohr

rprim   1.0000000000E+00  0.0000000000E+00  0.0000000000E+00
                       0.0000000000E+00  1.0000000000E+00  0.0000000000E+00
                       0.0000000000E+00  0.0000000000E+00  1.0000000000E+00
       
natom    14  
ntypat   3  
typat    1 1 1 1 2 2 3 3 3 3 3 3 3 3 
znucl    38 22 8  
     
xred     -3.4744227802E-18  5.0000000000E-01  9.7732944288E-02
                      -3.4744227802E-18  5.0000000000E-01  2.7726705571E-01
                       5.0000000000E-01 -3.4710334963E-18  9.7732944288E-02
                       5.0000000000E-01 -3.4710334963E-18  2.7726705571E-01
                      -3.4744227802E-18 -3.4710334963E-18  1.8750000000E-01
                       5.0000000000E-01  5.0000000000E-01  1.8750000000E-01
                      -3.4744227802E-18 -3.4710334963E-18  9.3765936754E-02
                      -3.4744227802E-18 -3.4710334963E-18  2.8123406325E-01
                       5.0000000000E-01  5.0000000000E-01  2.8123406325E-01
                       5.0000000000E-01  5.0000000000E-01  9.3765936754E-02
                       2.9796231887E-01  2.0203768113E-01  1.8750000000E-01
                       7.0203768113E-01  7.9796231887E-01  1.8750000000E-01
                       2.0203768113E-01  7.0203768113E-01  1.8750000000E-01
                       7.9796231887E-01  2.9796231887E-01  1.8750000000E-01

pp_dirpath "~/pseudo/PBEsol/"
pseudos "Sr.psp8, Ti.psp8, O.psp8"

# Basis Set
ecut    36.0 
ecutsm    0.5     

nband   70
nbdbuf 14

ngkpt   16 16 1
nshiftk 1
shiftk  0.0 0.0 0.0

# SCF Procedure
nstep   500       
diemac  6.0

I was using the cells that were fully converged to begin with (tolvrs 1.0d-10) in a fixed lattice strained state (optcell 0).

I am attaching my relaxation and phonon input/ouptut files here for reference.

Thanks
Dominic

phonon_qqq.abo (192.4 KB)
phonon_qqq.abi (2.4 KB)
latt_opt.abo (144.1 KB)
latt_opt.abi (1.9 KB)