Non-constant Vacuum potential

Question

I relaxed Na₂S₂ molecule in QuantumEspresso and plotted the planar average of potential energy across z axis so as to find the Vacuum potential. This is the output I got : As you can see the vacuum potential is not constant making it impossible to calculate it. While following the same procedure for Na₂S I got the following: as you can see the vacuum potential is a constant value. What could be a reason for what happened in the case of Na₂S₂? How can I solve it ?

Thank you!

The input file for the geometry optimization of Na₂S₂:

&CONTROL
calculation='relax'
title='graphene'
prefix='graphene'
restart_mode='from_scratch'
outdir='../tmp'
wf_collect=.true.
pseudo_dir='~/work/pot/'
/
&SYSTEM
ibrav = 4,
a = 14.58 ,
c= 20,
nat = 4
ntyp = 2
vdw_corr = 'dft-d3'
ecutwfc = 50.0 ,
ecutrho = 250.0 ,
occupations='smearing'
smearing='gaussian'
degauss=0.001
nbnd=30
/
&ELECTRONS
electron_maxstep = 200,
conv_thr = 1.0d-10 ,
/
&IONS
ion_dynamics = 'bfgs'
/

ATOMIC_SPECIES
Na 22.989 Na.pbe-spnl-kjpaw_psl.1.0.0.UPF
S 32.065 S.pbe-n-kjpaw_psl.1.0.0.UPF

ATOMIC_POSITIONS (angstrom)
Na           -2.1578688813       -0.0711007913        0.0959862818
S             0.0279507567        0.7976264059        1.1138133631
S            -1.4993815540        2.3462755509        0.6717950344
Na            0.8998813708        2.9052562024       -0.0587466613


 K_POINTS automatic
2 2 1 0 0 0

Answer 1

As mentioned by Phil Hasnip in the comments, placing a $\mathrm{Na_2S_2}$ molecule under periodic boundary conditions is akin to simulating a 3d array of these molecules. If the molecule has a nonzero dipole moment, this give rise to an electric (dipole) field between adjacent 2d "planes" of molecules. Unless you take care to orient the dipole of your molecule along a specific cartesian axis, it can point in any direction of the cell. I.e. the effect you are seeing along the z direction may as well exist along x and y with different magnitudes proportional to the component of the molecular dipole in that direction.

As with a parallel plate capacitor, at sufficient distance from the 2d "planes" of molecules the electric field is roughly constant, giving rise to the linear region of the potential in your plot. The purpose of the dipole correction is to place a counter-dipole of opposite magnitude at a position sufficiently far away from the molecule in order to "flatten" the potential curve. I'm not quite sure what you mean by "when I add dipole correction gave me two vacuum potentials". When done right (see e.g. this example by Christoph Wolf), you should simply obtain a flat potential.

By the way, if you are just interested in the value of the vacuum potential, you could place the molecule in the center of the cell and simply take the value of the potential in a corner of the simulation box. This neglects the polarizing effect of the electric field on the electron cloud in the molecule but should converge to the vacuum potential with increasing cell size.

Finally, you don't state the purpose of this calculation (and there may be valid reasons for you to use Quantum ESPRESSO) but if it is to compute quantities like the ionization potential or the electron affinity of a molecule, don't use a code with periodic boundary conditions. This is what quantum chemistry codes are for, and there are both great commercial and free implementations.