Lack of translational symmetry using using Be_sv with PBE

[svb]

I am getting unphysical results when using Be_sv PPs with PBEsol. Suspect a serious underlying bug. Seen using vasp 5.3.5, please alert if this relates to a bug fixed in a later release.

==Summary===

Looks like it boils down to a lack of translational symmetry: for two identical geometries (rigid shift in the direction of a first unit vector), I get 2.3meV difference despite using a very high precision. This is the first example below.

The second example below is a complete lack of acoustic modes in the phonon spectra. While in LDA and PBE the lowest 6 modes involved a triply degenerate acoustic mode and a triply degenerate soft mode, with the eigenvectors of the acoustic modes proportional to sqrt(m) and the eigenvectors of the soft mode ~sqrt(1/m) to a reasonable accuracy [whether the soft mode was real or imaginary, depending on the volume], in PBEsol with Be_sv PP there appeared a strange mode with all atoms displacing in the same direction (as in the acoustic mode), but the eigenvector significantly deviating from ~sqrt(m), and the imaginary frequency as high as 1.4THz (still being the lowest absolute value frequency)! Clearly this is a mixture of an acoustic and the soft mode, but those modes should not ever mix (even if degenerate) as long as there is a translational symmetry!

Below are the details of the setup and the unphysical results for those two cases.

======================================================================
==================---- First example ---==================================
Relaxation using pre-relaxed PBEsol geometry is performed for two equivalent poscars:
First POSCAR:
-------------------------------------------------------------------
Rocksalt BeO cube
1.00000000000000
3.6229349637732597 -0.0000000000000000 0.0000000000000000
-0.0000000000000000 3.6229349637732597 -0.0000000000000000
0.0000000000000000 -0.0000000000000000 3.6229349637732597
Be O
4 4
Direct
-0.0000000000000000 0.0000000000000000 -0.0000000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 -0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 -0.0000000000000000
-0.0000000000000000 -0.0000000000000000 0.5000000000000000
---------------------------------------------
Second POSCAR: same except displaced by 0.01 fractional in x direction
-------------------------------------------------------------------
Rocksalt BeO cube
1.00000000000000
3.6229349637732597 -0.0000000000000000 0.0000000000000000
-0.0000000000000000 3.6229349637732597 -0.0000000000000000
0.0000000000000000 -0.0000000000000000 3.6229349637732597
Be O
4 4
Direct
0.0100000000000000 0.0000000000000000 -0.0000000000000000
0.0100000000000000 0.5000000000000000 0.5000000000000000
0.5100000000000000 0.0000000000000000 0.5000000000000000
0.5100000000000000 0.5000000000000000 -0.0000000000000000
0.5100000000000000 0.5000000000000000 0.5000000000000000
0.5100000000000000 0.0000000000000000 0.0000000000000000
0.0100000000000000 0.5000000000000000 -0.0000000000000000
0.0100000000000000 -0.0000000000000000 0.5000000000000000
--------------------------
KPOINTS:
---------------------------
KPOINTS file
0
Gamma-centered
8 8 8
0 0 0
-------------------------
INCAR:
--------------------------
ISPIN = 1
LSORBIT = .FALSE.
PREC = High
ALGO = Normal # Normal (Davidson)
EDIFF = 1e-6
ISMEAR = -1
SIGMA = 0.200000 # eV
IBRION = 2
NSW = 50
POTIM = 0.2500000
ISIF = 3
LCHARG = .FALSE.
LWAVE = .FALSE.
LELF = .FALSE.
LVTOT = .FALSE.
LVHAR = .FALSE.

ADDGRID = .TRUE.
NPAR =1
GGA=PS #PBEsol
----------------------------
POTCAR versions:
-----------------
TITEL = PAW_PBE Be_sv 06Sep2000
TITEL = PAW_PBE O 08Apr2002
-------------------
Note using PREC=High, ADDGRID, and 8x mesh for 8-atom cell (for such large-gap insulators, already 4x mesh is usually sufficient to get well-converged energy differences and phonon frequencies). Nevertheless, the energies of the two equivalent structures are 2.3meV different:
grep F OSZICAR1
----------------------
1 F= -.55868584E+02 E0= -.55868584E+02 d E =-.558686E+02
2 F= -.55868604E+02 E0= -.55868604E+02 d E =-.201559E-04
3 F= -.55868564E+02 E0= -.55868564E+02 d E =0.191161E-04
4 F= -.55868604E+02 E0= -.55868604E+02 d E =-.201376E-04
----------------------
grep F OSZICAR2
----------------------
1 F= -.55866279E+02 E0= -.55866279E+02 d E =-.558663E+02
2 F= -.55866296E+02 E0= -.55866296E+02 d E =-.173921E-04
3 F= -.55866263E+02 E0= -.55866263E+02 d E =0.162595E-04
4 F= -.55866298E+02 E0= -.55866298E+02 d E =-.188656E-04
--------------------
Why is the translational symmetry violated?

======================================================================

==================---- Second example ---==================================

Phonon spectra: I tried calculating varying IBRION, ISYM, and used several cell volumes. Each time, as long as I use Be_sv with PBEsol, the mode that is closest to being acoustic clearly violates the translational symmetry [eigenvectors do not scale with sqrt(m)], and the smallest-absolute-value frequency is very different from zero. One particular example is given below.
POSCAR
-----------
Rocksalt BeO cube
1.00000000000000
3.6486000000000001 0.0000000000000000 0.0000000000000000
0.0000000000000000 3.6486000000000001 0.0000000000000000
0.0000000000000000 0.0000000000000000 3.6486000000000001
Be O
4 4
Direct
0.0000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.5000000000000000
0.5000000000000000 0.5000000000000000 0.0000000000000000
0.5000000000000000 0.5000000000000000 0.5000000000000000
0.5000000000000000 0.0000000000000000 0.0000000000000000
0.0000000000000000 0.5000000000000000 0.0000000000000000
0.0000000000000000 0.0000000000000000 0.5000000000000000
--------------------------
INCAR
-----------------
PREC = Accurate
ENCUT = 500
IBRION = 8
EDIFF = 1.0e-08
IALGO = 38
ISMEAR = -5
LREAL = .FALSE.
ADDGRID = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
LEPSILON = .TRUE.
GGA=PS #PBEsol
------------------------------
KPOINTS
-----------------
KPOINTS file
0
Gamma-centered
8 8 8
0 0 0
------------------------
POTCARs
-------------------
TITEL = PAW_PBE Be_sv 06Sep2000
TITEL = PAW_PBE O 08Apr2002
-------------------

Results:
there is no zero-frequency mode:
18 f = 9.248672 THz 58.111123 2PiTHz 308.502496 cm-1 38.249448 meV
19 f = 4.319060 THz 27.137456 2PiTHz 144.068338 cm-1 17.862204 meV
20 f = 4.319060 THz 27.137456 2PiTHz 144.068338 cm-1 17.862204 meV
21 f = 4.319060 THz 27.137456 2PiTHz 144.068338 cm-1 17.862204 meV
22 f/i= 1.393313 THz 8.754441 2PiTHz 46.475904 cm-1 5.762280 meV
23 f/i= 1.393313 THz 8.754441 2PiTHz 46.475904 cm-1 5.762280 meV
24 f/i= 1.393313 THz 8.754441 2PiTHz 46.475904 cm-1 5.762280 meV

Moreover, the eigenvector of the lowest mode (i*1.4THz) is NOT PHYSICALLY ALLOWED at a Gamma-point:
X Y Z dx dy dz
0.000000 0.000000 0.000000 -0.015525 -0.161782 -0.010574
0.000000 1.824300 1.824300 -0.015525 -0.161782 -0.010574
1.824300 0.000000 1.824300 -0.015525 -0.161782 -0.010574
1.824300 1.824300 0.000000 -0.015525 -0.161782 -0.010574
1.824300 1.824300 1.824300 -0.045060 -0.469576 -0.030692
1.824300 0.000000 0.000000 -0.045060 -0.469576 -0.030692
0.000000 1.824300 0.000000 -0.045060 -0.469576 -0.030692
0.000000 0.000000 1.824300 -0.045060 -0.469576 -0.030692
Note the Be and O eigenvectors differ by nearly a factor of 3, while the ratio should only be sqrt(m_O/m_Be)=1.31


For a comparison, using PBE instead of PBEsol and regular "PAW_PBE Be" instead of Be_sv, the ratio of the acoustic eigenvector components is reasonably close (though still far from being equal – a very annoying persistent problem in VASP!) to the expect ratio:
X Y Z dx dy dz
0.000000 0.000000 0.000000 -0.203310 0.210903 -0.056909
0.000000 1.824300 1.824300 -0.203310 0.210903 -0.056909
1.824300 0.000000 1.824300 -0.203310 0.210903 -0.056909
1.824300 1.824300 0.000000 -0.203310 0.210903 -0.056909
1.824300 1.824300 1.824300 -0.273321 0.283529 -0.076506
1.824300 0.000000 0.000000 -0.273321 0.283529 -0.076506
0.000000 1.824300 0.000000 -0.273321 0.283529 -0.076506
0.000000 0.000000 1.824300 -0.273321 0.283529 -0.076506

Similar results (except very different frequencies) using LDA and "Be". No clear signature of a problem using PBEsol and "Be", albeit the frequencies unexpectedly exceed those from both LDA and PBE, so not clear if those runs are problem-free.

What causes the unphysical mode when using Be_sv+PBEsol? (and while we are on it, would be great to learn what causes the typically ~3% deviation of the phonon eigenvectors from the symmetry-required ratios in virtually all VASP phonon calculations for all materials?)

Thank you very much!

[HechtSin]

Hi,

I think I might meet similar problem when I use "Cs_Sv" + PBEsol. The acoustic modes at zone-centers are not zeros.
Hope they can fix this.

[martin.schlipf]

While the real crystal is translational invariant, the numeric depiction of it is not. Because we use a finite real space grid, the relative position of the atoms to these grid points changes when you move all the atoms by a certain amount. Of course this can be reduced by increasing the energy cutoff ENCUT and thereby using a denser grid, but a small nonzero value for the acoustic phonons is expected.