# Density functional description of long-range electron Coulomb interactions
in bulk SnS

Stefanos Giaremis  Joseph Kioseoglou  Eleni Chatzikyriakou<EMAIL_ADDRESS>School of Physics, Aristotle University of Thessaloniki, 54124 Thessaloniki,
Greece

###### Abstract

A high-throughput benchmarking technique for testing the performance of
different exchange-correlation functionals and pseudopotentials is proposed
and applied to bulk SnS. It is shown that, contrary to the popular view that
the local density approximation can best describe layered materials, a
semilocal pseudopotential with a functional having a gradient dependence
better described lattice vectors and ‘tetragonicity’ of the lattice. We
classify the pseudopotentials based on this value and show that the
participation ratio of maximally localized Wannier functions follows the
theory which states that more distorted structures have higher anti-bonding
hybridization as stabilizing factor. In order to classify pseudopotentials,
the local and nonlocal potential contributions to the dynamical Born effective
charges are taken for each pseudopotential. Finally, a strategy is proposed
for learning exchange-correlation functionals based on the distinction between
short and long range parts of the Kohn-Sham potential.

##  I Introduction

Significant advances in first principles studies have recently opened new
opportunities for the examination of quantum phenomena in physics and
chemistry at a very detailed level and reduced cost [1, 2]. One important
issue that arises is the assessment of each technique, in terms of various
factors such as speed and accuracy [3].

High-throughput methods is another category that has gained popularity [4, 5,
6]. Beside analytical solutions used to screen large sets of proposed
compounds [7], other works use the same pseudopotentials, or pseudopotential
families, applied to all examined compounds, in order to find the desired
solutions. Another alternative is to use high-throughput in conjunction with
machine learning techniques, with different applications including the
detection of ferroelectric materials and testing of appropriate descriptors
[8]. Here we follow a different approach, in which we use high-throughput
techniques in order to benchmark the pseudopotentials themselves. We do this
following a physically-motivated procedure that lays the grounds for machine
learning dynamical lattice properties, without the use of separate classical
molecular dynamics run, but in combination with first principles electronic
properties in Density Functional Theory (DFT). This combined calculation is
also more convenient in the case of layered compounds, where the long-range
interactions play a significant role both for both structure and properties,
such as in twisted multi-layers.

Due to limitations of the exchange-correlation (XC) functionals in
realistically simulating interactions of different range orders [9], internal
properties such as binding energies, polarization and phonon effects can
present a qualitative difference in the amount of error introduced compared to
more ‘combined’ quantities, such as, lattice constants. In SnS, the internal
pressure and forces at different directions that result from the involvement
of the lone pairs, makes it hard to give accurate quantitative results. For
example, the free-standing 2D systems, can be close to a phase transition and
therefore predicted to be stabilized with the addition of a substrate, but it
is uncertain whether a different pseudopotential (PP) would produce similar
results [10]. In this work, we benchmark pseudopotentials for their ability to
describe this process in SnS, however, our results can be generalized to other
types of compounds containing lone pairs, such as post-transitions metal (Tl,
Pb, B), or other types of chalcogenide-based compounds (i.e., GeSe, SnSe).

###  I.1 Crystal symmetry

SnS has many intriguing properties in its 2D form, however, the bulk form has
also shown important properties that include room-temperature valley physics
[11, 12, 13]. In this work, we examine the bulk form in order to avoid the
depolarization field effects that arise from the periodic nature of Density
Functional Theory. Although there is enough degree of general consensus over
its crystal structure, the details have been controversial until very recently
[14, 15].

Sn(II) compounds like SnS form covalent bonds with mixed ionic character [16].
In this text, we follow the notation of the direction being perpendicular to
the separate layers, and the x and y being the in-plane directions (fig. 1).
In this configuration, the bonds are formed between Sn and S orbitals of the
same layer in both the in-plane and out-of-plane directions. Sn lone pairs
extend in the void between, what appears to be, weakly bonded layers.

Figure 1: Crystallographic direction notation used in the text.

It was realized from early on that lone pair formation was important for the
lattice structure. In the revised lone pair (RLP) model [17, 18], this view
persists[19]. Distortions of the lattice are due to the hybridization of the
Sn(5p) with the Sn(5s)+S(3p) anti-bonding orbitals (Figure 2). It is generally
assumed that the more the energies of the latter overlap, the more strongly
they interact, producing stronger anti-bonding states, which hybridize with
Sn(5p) states. This creates a higher distortion, and a larger interlayer
distance. The resulting lattice is less symmetrical, with angles away from 90o
(e.g. SnO). Weaker interactions (SnS, SnSe), and eventually breaking of such
bonds (SnTe), leads to more symmetrical structures [20].

Figure 2: Energy level diagram for Sn(II) bonding according to the RLP model
[17, 18]. The Sn 5p states hybridize with the (Sn 5s - anion p) anti-bonding
states.

Similar symmetry-breaking transitions also occur as a function of temperature
of bulk SnS from the higher-temperature, higher-symmetry Cmcm spacegroup to
the the low-temperature, lower-symmetry Pnma [21, 22], and the same orbital
hybridization model is believed to be behind this mechanism [23].

While these facts are well-known, and many authors choose to not include vdW
effects in their calculations [24], what is generally not mentioned often is
the important role the stereochemically active lone pairs play at stabilizing
the structure. These anti-bonding orbitals, produce a stabilizing force that
is stronger than the weak vdW forces [25, 26], and they are also the reason
single-crystal 2D SnS is hard to achieve experimentally on a large scale [27,
28, 29, 30, 31, 32].

It is also worth mentioning that recent studies showed that the RLP model is
not fully applicable to the case of SnTe, as well as other cases such as
(PbSe, PbTe, GeTe). In these cases, the type of bonding has been named
metavalent, and is intermediate between metallic and covalent or ionic bonding
[33]. While, in principle there materials should also have stronger interlayer
bonds than vdW type, the atoms share approximately one electron between them
[34]. The region where the electron extends and its localization properties
depend on many intrinsic (e.g. dimensionality) as well as extrinsic factors.

##  II Density functional approximation

The starting point of our analysis is the Kohn–Sham potential, whose general
expression is,

$V_{KS}=V_{ion}+V_{H}+V_{XC}$ (1)

where the terms correspond to the electron-ion interaction, Hartree and XC
potential energies respectively. This, along with the Hohenberg-Kohn theorem,
relates the energy with the density of the system and allows for the
separation of the contributions of each type of potential. The kinetic energy
of the electrons in inferred from the charge density that is produced from
this potential energy after the self-consistent procedure.

###  II.1 Exchange-Correlation functional

Starting with VXC, we note the important difference between the performance of
the Generalized Gradient Approximation (GGA) functional and the Local Density
Approximation (LDA). Short-range electron correlation effects are explicitly
included in the DFT functionals, however, long-range correlations are also
thought to be included through the self-interaction error. The underestimation
of the approximated exchange hole (that depends on the exchange functional) is
equivalent to a re-construction of the bands around the reference electron and
mimicking nondynamic, long-range correlation effects [9].

In general, the LDA exchange hole is diffuse and overestimates the correlation
energy. It produces larger self-interaction errors than GGA for the core and
lone pair orbitals. At the same time, the GGA hole is more localized than the
HF exchange hole in the bond and lone pair regions due to the Coulomb self-
repulsion being over-compensated [9].

It is generally established that the majority of the GGA functionals
exaggerate bond length, while LDA ‘overbinds’. LDA functionals are thought to
better describe layered materials, that is, a description that involves
regions of low electron density [35], and this has been adopted for SnS as
well [22, 36, 37]. There is a second factor that affects accuracy for
materials with polarization properties, that is the insensitivity of the local
functional to the polarization, when the exact XC functional is polarization
dependent [38, 39, 40]. This, however, becomes a problem only when the
electronic response becomes equally significant as the ionic response, which
is many times not the case in ferroelectric oxides.

As described in the previous section, stabilization of SnS is a combination of
the hybridization between the bonding and the anti-bonding orbitals that
reside at distances where the approximation methods produce qualitatively
different results. This means that there exist contributions from both bonding
and non-bonding regions, while the long-range vdW interaction terms play a
secondary role. At the same time, it pays to consider the lattice and
electronic contributions to the dielectric response of the system, and the
quantities produced by each approximation method.

In the self-interaction-corrected (SIC) approximation, the exchange-
correlation energy of a single, fully occupied orbital is taken to cancel
exactly the self-Coulomb energy (the functional of the electron number
density) [41]. While this has been shown to improve some total and ionization
energies in some systems, it has also been shown to worsen the results in a
plethora of systems that need the self-interaction error in order to
compensate for the lack of long-range force description in DFT [9]. In all
cases, the SIC effects are reduced as a function of distance from the ion
[41]. In our case, for the 5s and 5p orbitals of Sn this effect might become
negligible.

In order to keep the computational time and general simulation complexity low,
only approximations that belong to the first rungs of the Jacob’s ladder were
used. For LDA, Slater-type exchange [42] and Perdew–Wang (PW) correlation
energy [43] are assumed, while for GGA, the Perdew, Burke, Ernzerhof (PBE)
gradient correction to the exchange and correlation [44]. The solids-oriented
PBEsol, with the diminished gradient dependence [45] was also considered, as
well as the SIC Perdew Zunger (PZ) LDA functional [41].

###  II.2 Electron-ion potential

Four basic types of approximations of the electron-ion interaction potential
energy, Vion, have been examined. Starting with the semilocal (SL) PP,

$\begin{split}V_{ion}=V_{L}+\sum_{l}V_{SL,l}\hat{P}_{l}\\\
V_{SL,l}=V_{l,NL}-V_{L}\end{split}$ (2)

where the subscript L denotes local pseupotential, NL the non-local potential
and l the angular momentum component of the wave-function. For a nonlocal
ionic potential, a different pseudopotential, Vl, should be acting on each
angular momentum component,

$V_{NL}=\sum_{lm}\bra{\Psi_{lm}}V_{l}\ket{\Psi_{lm}}$ (3)

where $\Psi_{lm}$ are the spherical harmonics. Semilocal PPs essentially
become non-local only beyond a cutoff radius, a process which allowed the
plane-wave basis set to be expanded to the nonlocal description by reducing
the computational cost, but is also a source of error.

SL PPs are norm-conserving, that is, the norm of the pseudo wavefunction is
conserved inside the core region [46]. Another form of the norm-conserving PPs
was given by Rappe, Rabe, Kaxiras, and Joannopoulos (RRKJ) [47, 48] who
expanded each pseudo-wave-function as a linear combination of basis functions.
The implementation of this PP in the present work is ultrasoft (US).

For PAW PPs, we use the KJ derivations [49]. In the PAW formalism, the atomic-
based projectors $\beta$ participate in the KS potential during the self-
consistent cycle through
$\sum_{i}\braket{\psi_{i}}{\beta_{n}}\braket{\beta_{m}}{\psi_{i}}$, with
$\psi_{i}$ the calculated wavefuction [50].

In Ref. [51], Hamman also proposed the optimized version of the norm-
conserving Vanderbilt PPs, which rely on the multiple projector scheme
proposed by Vanderbilt for passing from norm-conserving to ultrasoft PPs.
These are the Optimized Norm Conserving (ONCV) used here. The details of the
implementations are not important for this work, because their differences can
be captured in the local and non-local definition of the Born Effective
charges given in the next section.

###  II.3 Range of interactions

The existence of ferroelectricity in displacive ferroelectrics, results from
the interplay between short-range and long-range interactions [52, 53]. That
is, only short range interactions are required to support the high-symmetry
paraelectric phase, whereas long-range interactions are required, in order for
the ferroelectric phase to persist in any configuration (when raising the
temperature, thinning down to a few layers or doping). If the long-range
interactions are somehow screened, ferroelectricity can, in principle, be
suppressed, even above the transition temperature [52]. Here we examine the
effects of pseudopotential in the description of lattice dynamical properties,
as they arise from differences in the short and long range interactions of the
electrons and ions of SnS.

The ferroelectricity results from the broken centrosymmetry in odd number of
2D SnS. In the literature, 2D SnS is considered displacive [54], while the
defining factor of the ferroelectric/paraelectric properties is the short-
range interaction between the lone pairs of electrons [22]. Here we review the
quantities involved in the calculations in a first principles set-up.

Taking each term in equation 1: VXC is short ranged, VH has a short and a
long-range part and Vion is long (infinite) range. The later is further
divided in local and non-local part, VL and VNL respectively. Of the two, the
non-local potential is short ranged, while the local potential has both short-
range (SR) and long-range(LR) contributions [55, 56]. Therefore,

$\begin{split}V_{KS}^{SR}=V_{L}^{SR}+V_{H}^{SR}+V_{XC}+V_{NL}\\\
V_{KS}^{LR}=V_{L}^{LR}+V_{H}^{LR}\end{split}$ (4)

In DFT codes, the convention is that the SR and LR parts are dictated by the
cut-off radius during construction of the PP, and then participate in the
Ewald summation with the former being calculated in real space, while the
later in reciprocal space.

###  II.4 Dynamical properties

Independent of the method used to derive them, the Born effective charges,

$Z^{*}=-\partial\mathbf{F}/\partial\mathscr{E}$ (5)

can be divided so as to only take into account contributions of a specific
potential. The forces are related to Vion,

$F_{\alpha,i}\propto\int_{\Omega}n\left(\mathbf{r}\right)\frac{\partial
V\textsubscript{ion}}{\partial\mathbf{u}_{\alpha,i}}$ (6)

for an atom $\alpha$ in the direction i and $\Omega$ the volume of the unit
cell. It should be noted that there is also an ion-ion contribution, as well
as a kinetic energy contribution in the calculation of the energy of the
system. The kinetic energy is short-ranged, but it is not included as a
separate term in the Kohn-Sham potential, but rather deduced from other
quantities during the self-consistent procedure. The ion-ion contribution only
depends on the atomic numbers of the elements and therefore should be similar
for all PPs.

##  III Relaxing forces

We performed the calculation in the Pnma phase of bulk SnS, with 4 Sn and 4 S
atoms. The zigzag and armchair directions are along the a and b axis,
respectively. First we calculated the relaxed bulk structures using each PP.
The results are presented in Table 1. They are listed in increasing value of
the a/b ratio. The value that most closely matches the experiment was given by
the SL PBE. The rest of the approximations gave values either lower or higher,
but with a significant difference from the PBE SL.

The LDA PPs produced higher ratios than GGA. This can be explained immediately
by the property of LDA to overbind. In fact, the more tetragonal structures
also have smaller in-plane lattice vectors, with the armchair direction
showing more pronounced reduction. The values given by the reduced gradient
dependence PBEsol, lie between those of LDA and PBE. Interestingly, the self-
interaction corrected LDA (PZ) produced identical a/b ratios with pure LDA
(PW), however, there was one significant difference: PW reduced the inter-
atomic distances of non-bonded Sn-S pairs, reaching bonding lengths, which
effectively changed the coordination of the atoms to five-fold. The only other
difference is an accompanying slight reduction in the $c$ lattice constant
produced with the PW. Overall, the SIC, had the effect of slightly elongating
the interatomic distances of the structure, compared to the nonphysical result
produced by the LDA.

The armchair (b) direction has the largest deviation between pseudopotentials,
which can be intuitively understood from the absence of bonds in this
direction that results in lower electron density, in which DFT has been
notoriously bad at approximating.

The scenery changes slightly when we consider the out-of-plane lattice
vectors. Here, the SL PP performed only marginally better than LDA, with the
PBEsol giving the best estimate. Again the distances between atoms (interlayer
this time) produced by LDA, are smaller than those of PBE. This result shows
that, contrary to popular view, the interlayer distances given by LDA are not
necessarily better than those of GGA. The PBE results of SL and PAW for the c
lattice vectors improved with the addition of vdW interaction, however, they
remained qualitatively similar.

Table 1: Primitive lattice constants along the a, b and c directions after relaxation of bulk SnS without vdW forces included, compared to experimental measurements. The value that most closely matched the experiment is given in blue. The results are arranged in increasing a/b ratio.  |  a [Å] (%) |  b [Å] (%) |  c [Å] (%) |  a/b (%)
---|---|---|---|---
exp. [27] | 3.987 | 4.334 | 11.199 | 0.919
PBE PAW | 4.027(+1.0) | 4.454(+2.7) | 11.464(+2.3) | 0.904(-1.7)
PBE ONCV | 4.029(+1.0) | 4.450(+2.6) | 11.426(+2.0) | 0.905(-1.5)
PBE US | 4.028(+1.0) | 4.440(+2.4) | 11.441(+2.1) | 0.907(-1.3)
PBE SL | 4.023(+0.91) | 4.360(+0.5) | 11.362(+1.4) | 0.922(+0.3)
PBEsol US | 3.989(+0.05) | 4.238(-2.1) | 11.107(-0.8) | 0.941(+2.2)
PBEsol PAW | 3.989(+0.07) | 4.234(-2.2) | 11.106(-0.8) | 0.942(+2.4)
PZ US | 3.961(-0.6) | 4.189(-3.3) | 10.995(-1.8) | 0.945(+2.7)
PZ PAW | 3.961(-0.6) | 4.187(-3.3) | 10.999(-1.7) | 0.945(+2.7)
PW ONCV | 3.957(-0.7) | 4.182(-3.4) | 10.975(-1.9) | 0.946(+2.8)

The case of under-convergence of wavefunction kinetic energy cut-off,
$e_{\mathrm{T}}$, in terms of the resulting internal pressure was also
examined. For some ferroelectrics, transitions have been reported at -5 GPa
[57]. For SnSe, it has been proposed that the Curie temperature can be tuned
up to a few hundred Kelvin, only by the application of a small strain, on the
order of 1% deviation of the lattice constant in the armchair direction [54].
The differences in different pseudopotentials in describing internal
properties such as pressure have been widely examined in the literature for
different types of ferroelectrics. For example, in PbTiO3, a supertetragonal
structure was obtained either with GGA PPs, or with LDA at negative pressure
[58].

Figure 3: Internal pressure vs. a/b lattice constant ratio (top) and c lattice
constant (bottom). Different points represent different kinetic energy cutoff,
$e_{\mathrm{T}}$ for the wavefunctions. The energy cut-off for charge density,
$e_{\mathrm{\uprho}}$, is 4 $\cdot\ e_{\mathrm{T}}$ in all cases, except for
ultrasoft pseudopotentials, where $e_{\mathrm{\uprho}}=8\cdot e_{\mathrm{T}}$.
The solid vertical line shows the experimental value. Unless stated as ‘PAW’,
the plane-wave method is used. ‘PW’ refers to the LDA PW functional [43].

The results of the relaxation are shown in fig. 3. Some approximations (i.e.
PBEsol) generally require higher values of $e_{\mathrm{T}}$. Lower values
produced higher negative internal pressures, but most converged to positive
values close to zero at higher $e_{\mathrm{T}}$. The values of a/b ratios, as
well as c lattice constants only slightly changed for the different values of
pressure, but the changes were significant for PBE PAW, which required
$e_{\mathrm{T}}\ >$ 75 Ry to achieve convergence. Some authors tend to use the
experimental lattice constants for their self-consistent calculations. In our
simulations, there were noticeable internal pressures when the lattice
constants used for the self-consistent cycles did not equal to those of the
relaxed structure, so we did not adopt this technique anywhere in this work.

Figure 3 (bottom), further clearly shows that lower a/b ratios produce higher
interlayer distance, as is intuitively expected for more distorted structures.

##  IV Effects of bonding on structure

###  IV.1 Projected Density of States

To analyze the effects of bonding on the final structure produced, we used
Density of States calculations, projected on selected orbitals and sums
thereof (PDOS). The results are shown in fig. 4. This type of representation
does not distinguish between bonding and anti-bonding states, however, it
reveals details on the overlap of certain orbitals in their energy range. We
present only the results of three representative PPs. For the rest, the
results are almost overlapping to each case presented here, and their
similarities follow the trend presented in fig. 3. They are all included in
the supplementary material.

We follow the notation of Ref. [20] and divide the valence energy range in
regions I, II and III (fig. 4). As it is also visualized in fig. 4(a), region
I is the bonding region, region III is the anti-bonding region. In region II,
there exist almost none of the Sn(s) states. The energy difference between
them and the Sn(p) states, does not allow the two to directly hybridize. They
interact through mediation from the S(p) states [20].

Figure 4: Orbital-projected density of states for (a) the p orbitals of S and
(b) the s orbital of Sn and (c) the p orbitals of Sn, using three
representative pseudopotentials. The dashed lines in (b) show the centroid of
region I for each case.

Theory states that stronger anti-bonding hybridization produces larger
distortions, and lower a/b ratios. We do not present any quantitative results
on the integrals here, as it is not easy to set clear boundaries between the
regions. We leave this for the next section, where we present hybridization
based on the Wannier results. Here, we focus first on the energy ranges.

In region I, the S(p) states have almost matching energy overlap for each PP
used, except for a slight spread-out in energy with PBEsol. At the same time,
the S(p) states for PBE PAW have higher density in region III, getting
progressively lower for PBE SL and PBEsol. This shows that there are less S(p)
states in the anti-bonding region as we go towards more symmetrical
structures.

By looking at the centroid of the Sn(s) states in Region I in comparison with
the peak at the top of region II S(p) states, we can observe the energy
distance between them, which points to the amount of hybridization [20]. The
energy difference is roughly constant, and increases slightly for the more
symmetrical structure produced by PBEsol PAW, revealing a less hybridized
character of the states. Overall, the density of S(p) in region II gets lower
as we go from more distorted (PBE PAW) to more symmetrical structures
(PBEsol), the same an in region III.

Contrary to what was observed for the S(p) states in regions III and II, the
Sn(s) and Sn(p) states in all regions are almost identical for the three PPs,
except for one significant difference: With PBE PAW, the Sn(s) states have a
higher density than with PBE SL in region I, while for PBEsol they spread to
lower energies.

To sum up, we observe that while there are more S(p) states in the anti-
bonding energy region with increasing distortion, Sn does not seem to follow
the same trend. There are, however, differences in the Sn(s) states in the
bonding region.

###  IV.2 Wannier functions

We convert the DFT states into maximally localized Wannier functions (MLWF),
starting from projected WF on the valence s and p orbitals of the Sn and S
atoms, resulting to overall 32 bands. For the case of PBE SL, all valence
states were included, for the rest of the PPs, the Sn 4d states were excluded
from the wannierization using a disentanglement window. The details of the
convergence and checking of the WF results are given in the methods section.

Another method to examine the state would be the Density Overlap Regions
Indicator [59], however, this did not produce valuable results, possibly
because it is mostly oriented towards non-bonding orbitals. A third method
exists, that is also suitable to high-throughput calculations, as it is devoid
of the need to specify an energy window and other parameters. It is based on
the ‘selected columns of the density matrix’ (SCDM) algorithm [60, 4] and is
able to derive Wannier functions at the highest electron density k-points.
Although such a method would allow for significant increase in automation, the
WFs produced depend on the localization properties of the states, that differ
for each PP, as shown below.

(a)

(b)

Figure 5: (a) Positions of the Sn and S atoms in the unit cell (b) Absolute
value of the hopping energy of intra-atomic bonding of Sn (left) and inter-
atomic bonding with the neighboring S atoms in the in-plane and out-of-plane
directions in the home unit cell (right). The integrals are taken from
projected WFs, before the minimization procedure. All Sn-S pairs shown are
bonded, except for Sn-S($\parallel$), which is a non-bonding pair.

In detail, the problem is related to the ‘projectability’ of the states.
Orbitals with more nodes, are harder to separate from the higher and lower
bands [4]. Therefore the lack of an energy window specification for this
method becomes a disadvantage. Indeed, we encountered both problems when
trying this method on SnS. To overcome this issue, either the parameters of
the SCDM, or the number of bands can be changed. We believe that the second
option would be more efficient as only a few conduction bands are required for
hybridization with the anti-bonding states, and the states become more
‘projectable’, the lower we are in energy (see supplementary information).

In order to automate the disentanglement procedure, the number of states
inside the energy window of interest was calculated from the first principles
bands calculation at each k point. The energy range was then re-adjusted so
that it becomes equal to the target WFs [6, 5]. Multiple checks were performed
regarding the validity and accuracy of the WFs. Without using the minimization
procedure, the bonding between the WFs, as revealed by the hopping integrals,
follows the trend in the bonding distance. In fig. 5(a) the position of the
atoms in the $y$ vs. $z$ directions after relaxation using DFT is plotted. The
Sn atoms are more displaced towards the direction of the lone pairs as we move
to more symmetric structures. The bond lengths expand and contract
accordingly.

In fig. 5(b) the hopping integrals between selected orbitals are plotted. The
atoms in the top and bottom layers are related to each other through inversion
symmetry. Similarly, for the atoms in each layer separately. The bonding
integrals are related to each other through such symmetries and the integrals
repeat accordingly, with the addition of a ($\pm$1 factor). Consequently, we
only show hoppings between orbitals of the top layer Sn and S atoms. It is
observed that as the symmetry of the structure increases, intra-atomic hopping
in the Sn atom decreases, in favor of inter-atomic hopping energy, which
increases. This reflects decrease in distances between the Sn and S atoms with
increasing symmetry.

Figure 6: Violin plots showing (top) the participation ratio of all the WFs in
all eigenenergies and (bottom) the participation ratio of the Sn(s,p) and S(p)
states in the antibonding region (III) including the first four conduction
band states. The two middle lines show the mean and median of the values.

To quantify that change in the localization properties of the states, and of
hybridization, the participation ratio of the MLWF is extracted. This is
defined as [61, 62],

$\mathrm{PR}\left(n\right)=\frac{\left[\sum_{i=1}^{N}|\braket{i}{n,\mathbf{k}=0}|^{2}\right]^{2}}{\sum_{i=1}^{N}|\braket{i}{n,\mathbf{k}=0}|^{4}}$
(7)

where N is the total number of MLWFs included in the wannierization and n is
the selected orbital state, while k = 0 denotes the home unit cell. The value
of the participation ratio decreases with increasing localization, and it more
properly reflects the number of orbitals over which hybridization of a state
occurs. In analogy to Ref. [63], we also define the ‘projected participation
ratio’, in the subspace Ws of the antibonding states. For this, the states of
eq. 7 were renormallized, so as to be defined only in the subspace $i\in
W_{s}$ of the 8 bands around the bandgap. The cut-off for this was taken so
that no states below region III in fig. 4 were taken.

The results are shown in fig. 6. As the symmetry increases, the states become
less hybridized overall (PR), or more localized. In the bottom panel (PPR),
only the Sn(s,p) and the S(p) states, which participate in the lone pair
model, are sampled. It is seen that the largest contribution to this trend in
localization comes from these states, projected in the anti-bonding subspace.
The RLP model is correctly recovered, which states that hybridization in the
distorted structures is higher. As was mentioned in section III, by using the
PW ONCV PP, the Sn-S distances reduce significantly, and the hybridization
between all neighboring atoms is wrongly emerging from the increase in
coordination. In this case, the high PR and PPR ratios reflect this
nonphysically higher coordination. The wannierization of this PP also did not
converge, producing high Im/Re ratios and higher spread values. The slow
convergence behavior of this PP, is another reflection of the high
hybridization level occurring. The rest of the PPs, showed very good to
excellent convergence behavior.

##  V Born effective charges

The dynamical atomic charges (Born Effective charges, Z*) are considered next
[64, 65]. The bulk periodic structure used here is not prone to depolarization
effects, so it is easier to compare between different pseudopotential
approximations. We calculated the Born effective charges of the Sn and S atoms
using the Berry-phase theory of polarization [66] in Quantum-Espresso [67,
68]. An electric field of 0.036 V/$\AA$ was used, which is the same order of
magnitude as the electronic gap in most cases ($\approx$0.01 V/$\AA$) [68].

Figure 7: Born effective charges of the Sn and S atoms separately, for each
PP. Sn1 and Sn2 refer to the bottom Sn atom of the top monolayer and the top
Sn atom of the bottom monolayer respectively. S2 is the in-plane neighbor of
Sn1 and S1 of Sn2. Figure 8: Components of the Born effective charges of the
Sn and S atoms and separating the local and non-local contributions of the
pseudopotential. Figure 9: Local potential for the Sn (top) and S (bottom)
atoms, as derived from the pseudopotential file.

The results are shown in fig. 7. We follow the same symmetry arguments as in
the case of the WFs and only present ‘non-repeated’ results. Our results for
ONCV-PBE are close to those presented for bulk SnS in Ref. [24] for the same
PP. The same trend is observed in the description of the Z* values, as in
Figure 3 for the ratios of the lattice constants to a very good approximation.
PBE PPs have lower values, which increase when passing to PBEsol and LDA
functionals. The two latter have almost equivalent Z* values.

The local (ZL) and non-local (ZNL) contributions for the atoms to the Z*
values were also extracted from the forces on the ions using equations 5 and
6. This is shown in fig. 8. In general, the magnitude of the contributions of
both VL and VNL in the Sn atoms is lower for the PBE than the rest of the
functionals, and higher for the S atoms. Especially in the case of SL, both
local and non-local contributions for Sn go to zero.

To be able to explain the differences in the functional, the information on VL
and VH around each atom, would suffice. This comes from the fact that the Z*
values are a result of the interplay of short-range and long range
interactions. Then, equation 4 indicates that the local and Hartree
contributions are enough to describe the long-range forces, if we are able to
divide the two contributions.

The visualization of VH is an involved process, and the results can not be
easily understood through the human eye, but would be more appropriate for
machine-learning methods. On the contrary, some insight can quickly be gained
directly from VL due to its rotational symmetry, and due to the fact that it
is available before the self-consistent procedure. By looking at its two parts
(short and long range) through the approximate division using the cut-off
radius, we can understand the differences at least between the PBE
functionals. This is shown in fig. 9, and in comparison to the LDA
functionals, in order to use them as a second reference for our study.

Notably, we observe that the VL for the Sn atom with the SL PP goes to zero
close to the core and is much shallower than the rest of the PPs. This
explains the $Z_{x,y,z}^{L}\approx 0$ observed in fig. 8. For PAW and US using
the PBE and PBEsol functionals, VL is equal. In these cases, it is the
functional that dictates the origin of the differences in the Z* values. The
same is true for the two PZ PPs examined, but only for the case of the S atom.
The last notable feature is in the case of PBE ONCV for the S atom. Here, the
large dip of VL is translated into large ZL values.

Overall, we were able to trace back many of the origins of the Z* values,
which makes our analysis interesting for physically-motivated derivation of
machine learning functionals and comparison with feature selection during
machine learning pseudopotentials for dynamical lattice parameters and
ferroelectric applications.

## Conclusions

We have successfully shown the dependence of the static and dynamic lattice
properties to different types of pseudopotentials.

The PPs could be categorized based on the a/b lattice constant ratio they
produced. The categorization also produced the same trend in the results for
the density of states, when approximately dividing the energies in regions
that belong to the bonding, intermediate and anti-bonding regions. This
reflects the dependence of the density functionals in the revised lone pair
model that dictates the stability of the layered structure to the lone pairs
in between the layers. The dependence was further revealed by the
participation ratio of the MLWFs, which point to higher hybridization for
higher distortion, and also showed that nonphysical bonding produced by the
density functional and PP combination had very bad convergence during
wannierization.

In terms of the dynamical lattice properties, the trend in the resulting Z*
values also followed the categorization that the PPs produced for the a/b
ratios. We traced back the properties of the local and non-local potential
contributions to the Z* values and also showed the connection of the former to
VL.

One very interesting direction of this is the use of the results to train a
model that learns dynamical lattice properties that can be measured
experimentally and then build the functional by using experimental results
[69, 70]. The combination of short-range interactions of the local potential
of the PP and the short-range Hartree interactions could then be used for this
purpose, aiming at displacive ferroelectric materials.

## Methods

DFT calculations were performed via the quantum espresso package [67]. The
pseudopotentials used were: 1) Ultra-soft, scalar-relativistic, Troullier-
Martins [46] PBE PAW with non-linear core correction [71] of Kresse-Joubert
type [49] (PBE PAW), 2) Ultra-soft, scalar-relativistic PBE by Rappe, Rabe,
Kaxiras and Joannopoulos (RRKJUS) [47, 48] with non-linear core correction, 3)
Ultra-soft, scalar-relativistic PBEsol RRKJUS with non-linear core correction
(PBEsol US) [45], 4) Ultrasoft, scalar-relativistic PBEsol PAW with non-linear
core correction (PBEsol PAW), 5) Semilocal, scalar-relativistic, Troullier-
Martins PBE, 6) Scalar-relativistic, optimized norm-conserving Vanderbilt PBE
pseudopotential (PBE ONCV), 7) LDA, scalar-relativistic, optimized norm-
conserving Vanderbilt PW [43] (PW ONCV), 8) Ultra-soft LDA, scalar-
relativistic PAW with non-linear core correction and self-interaction
correction (PZ PAW) [41] and 9) Ultra-soft, scalar-relativistic, LDA of
Troullier-Martins type (LDA US) with non-linear core correction [72].

The parametrization for all pseudopotentials is as performed in the PSlibrary
v1.0.0 by Dal Corso [72], except for the all ONCV type and the LDA US, which
were parametrized by Hamman [51, 73], while the PBE SL is taken from the
FHI98PP package (PBE SL) [74]. The forces acting on the atoms were relaxed to
a value below 2.5 meV/$\mathrm{\AA}$, however, the results with the chosen eT
did not change significantly when relaxing to a value up to three orders of
magnitude lower for most PPs. For the case were Van der Walls corrections were
included, this was done via the semi-empirical DFT-D3 method by Grimme [75].

The aiida software [4, 76] was used for deriving the projectabilities of the
states for the SCDM method. The convergence for the projected WFs (both with
and without localization) was set to be achieved when both the gauge invariant
and the gauge non-invariant parts of the total spread between successive
iterations was less than 10-9 \- 10-10 $\AA^{2}$, depending on the PP used.
The imaginary/real ratios were mostly zero or below 10-5, except for PZ PAW,
where convergence was not achieved with these settings. A k-point of 11x11x11
was used in all cases. TBmodels was used for the extraction of the hopping
integrals from the tight binding model of the wannier functions [6, 5].

##  VI Acknowledgment

Results presented in this work have been produced using the AUTH Compute
Infrastructure and Resources and was supported by computational time granted
from the Greek Research & Technology Network (GRNET) in the ‘ARIS’ National
HPC infrastructure under the project NOUS (pr012041).

## References

  * Hermann  _et al._ [2022] J. Hermann, J. Spencer, K. Choo, A. Mezzacapo, W. M. C. Foulkes, D. Pfau, G. Carleo, and F. Noé, Ab-initio quantum chemistry with neural-network wavefunctions,  (2022), arXiv:2208.12590 .
  * Carleo  _et al._ [2019] G. Carleo, K. Choo, D. Hofmann, J. E. Smith, T. Westerhout, F. Alet, E. J. Davis, S. Efthymiou, I. Glasser, S.-H. Lin, M. Mauri, G. Mazzola, C. B. Mendl, E. van Nieuwenburg, O. O’Reilly, H. Théveniaut, G. Torlai, F. Vicentini, and A. Wietek, NetKet: A machine learning toolkit for many-body quantum systems, SoftwareX 10, 100311 (2019).
  * Wu  _et al._ [2023] D. Wu, R. Rossi, F. Vicentini, N. Astrakhantsev, F. Becca, X. Cao, J. Carrasquilla, F. Ferrari, A. Georges, M. Hibat-Allah, M. Imada, A. M. Läuchli, G. Mazzola, A. Mezzacapo, A. Millis, J. R. Moreno, T. Neupert, Y. Nomura, J. Nys, O. Parcollet, R. Pohle, I. Romero, M. Schmid, J. M. Silvester, S. Sorella, L. F. Tocchio, L. Wang, S. R. White, A. Wietek, Q. Yang, Y. Yang, S. Zhang, and G. Carleo, Variational Benchmarks for Quantum Many-Body Problems, , 1 (2023), arXiv:2302.04919 .
  * Vitale  _et al._ [2020] V. Vitale, G. Pizzi, A. Marrazzo, J. R. Yates, N. Marzari, and A. A. Mostofi, Automated high-throughput Wannierisation, npj Computational Materials 6, 10.1038/s41524-020-0312-y (2020).
  * Gresch  _et al._ [2018] D. Gresch, Q. Wu, G. W. Winkler, R. Häuselmann, M. Troyer, and A. A. Soluyanov, Automated construction of symmetrized Wannier-like tight-binding models from ab initio calculations, Physical Review Materials 2, 103805 (2018), arXiv:1805.12148 .
  * Gresch  _et al._ [2017] D. Gresch, G. Autès, O. V. Yazyev, M. Troyer, D. Vanderbilt, B. A. Bernevig, and A. A. Soluyanov, Z2Pack: Numerical implementation of hybrid Wannier centers for identifying topological materials, Physical Review B 95, 1 (2017), arXiv:1610.08983 .
  * Olsen  _et al._ [2019] T. Olsen, E. Andersen, T. Okugawa, D. Torelli, T. Deilmann, and K. S. Thygesen, Discovering two-dimensional topological insulators from high-throughput computations, Physical Review Materials 3, 024005 (2019), arXiv:1812.06666 .
  * He  _et al._ [2021] J. He, J. Li, C. Liu, C. Wang, Y. Zhang, C. Wen, D. Xue, J. Cao, Y. Su, L. Qiao, and Y. Bai, Machine learning identified materials descriptors for ferroelectricity, Acta Materialia 209, 116815 (2021).
  * Polo  _et al._ [2003] V. Polo, J. Gräfenstein, E. Kraka, and D. Cremer, Long-range and short-range Coulomb correlation effects as simulated by Hartree-Fock, local density approximation, and generalized gradient approximation exchange functionals, Theoretical Chemistry Accounts 109, 22 (2003).
  * Malyi  _et al._ [2019] O. I. Malyi, K. V. Sopiha, and C. Persson, Energy, phonon, and dynamic stability criteria of two-dimensional materials, ACS Applied Materials & Interfaces 11, 24876 (2019), https://doi.org/10.1021/acsami.9b01261 .
  * Lin  _et al._ [2018] S. Lin, A. Carvalho, S. Yan, R. Li, S. Kim, A. Rodin, L. Carvalho, E. M. Chan, X. Wang, A. H. Castro Neto, and J. Yao, Accessing valley degree of freedom in bulk Tin(II) sulfide at room temperature, Nature Communications 9, 1455 (2018).
  * Chen  _et al._ [2018] C. Chen, X. Chen, Y. Shao, B. Deng, Q. Guo, C. Ma, and F. Xia, Valley-Selective Linear Dichroism in Layered Tin Sulfide, ACS Photonics 5, 3814 (2018).
  * Chatzikyriakou and Kioseoglou [2019] E. Chatzikyriakou and J. Kioseoglou, Emergence of valley selectivity in monolayer tin(ii) sulphide, Nanoscale Advances 1, 4863 (2019).
  * Burton and Walsh [2012] L. A. Burton and A. Walsh, Phase stability of the earth-abundant tin sulfides SnS, SnS2, and Sn2S3, Journal of Physical Chemistry C 116, 24262 (2012).
  * Biacchi  _et al._ [2013] A. J. Biacchi, D. D. Vaughn, and R. E. Schaak, Synthesis and crystallographic analysis of shape-controlled SnS nanocrystal photocatalysts: Evidence for a pseudotetragonal structural modification, Journal of the American Chemical Society 135, 11634 (2013).
  * Rundle and Olson [1964] R. E. Rundle and D. H. Olson, The Crystal Chemistry of Divalent Tin, Inorganic Chemistry 3, 596 (1964).
  * PrI [1967] _Progress in Inorganic Chemistry_ , Vol. 8 (1967) pp. 357–444.
  * Walsh  _et al._ [2011] A. Walsh, D. J. Payne, R. G. Egdell, and G. W. Watson, Stereochemistry of post-transition metal oxides: Revision of the classical lone pair model, Chemical Society Reviews 40, 4455 (2011).
  * Payne  _et al._ [2006] D. J. Payne, R. G. Egdell, A. Walsh, G. W. Watson, J. Guo, P.-A. Glans, T. Learmonth, and K. E. Smith, Electronic Origins of Structural Distortions in Post-Transition Metal Oxides: Experimental and Theoretical Evidence for a Revision of the Lone Pair Model, Physical Review Letters 96, 157403 (2006).
  * Walsh and Watson [2005] A. Walsh and G. W. Watson, Influence of the anion on lone pair formation in Sn(II) monochalcogenides: A DFT study, Journal of Physical Chemistry B 109, 18868 (2005).
  * von Schnering and Wiedemeier [1981] H. G. von Schnering and H. Wiedemeier, The high temperature structure of $\beta$-SnS and $\beta$-SnSe and the B16-to-B33 type A-transition path, Zeitschrift fur Kristallographie \- New Crystal Structures 156, 143 (1981).
  * Lebedev [2018] A. I. Lebedev, Ferroelectricity and piezoelectricity in monolayers and nanoplatelets of SnS, Journal of Applied Physics 124, 10.1063/1.5035419 (2018).
  * Li  _et al._ [2015] C. W. Li, J. Hong, A. F. May, D. Bansal, S. Chi, T. Hong, G. Ehlers, and O. Delaire, Orbitally driven giant phonon anharmonicity in SnSe, Nature Physics 11, 1063 (2015).
  * Dewandre  _et al._ [2019] A. Dewandre, M. J. Verstraete, N. Grobert, and Z. Zanolli, Spectroscopic properties of few-layer tin chalcogenides, JPhys Materials 2, 1 (2019).
  * Song and Lü [2018] H.-y. Song and J.-t. Lü, Density functional theory study of inter-layer coupling in bulk tin selenide, Chemical Physics Letters 695, 200 (2018).
  * Xu  _et al._ [2022] B. Xu, J. Deng, X. Ding, J. Sun, and J. Z. Liu, Van der Waals force-induced intralayer ferroelectric-to-antiferroelectric transition via interlayer sliding in bilayer group-IV monochalcogenides, npj Computational Materials 8, 10.1038/s41524-022-00724-8 (2022).
  * Ehm  _et al._ [2004] L. Ehm, K. Knorr, P. Dera, A. Krimmel, P. Bouvier, and M. Mezouar, Pressure-induced structural phase transition in the IV–VI semiconductor SnS, Journal of Physics: Condensed Matter 16, 3545 (2004).
  * Higashitarumizu  _et al._ [2020] N. Higashitarumizu, H. Kawamoto, C. J. Lee, B. H. Lin, F. H. Chu, I. Yonemori, T. Nishimura, K. Wakabayashi, W. H. Chang, and K. Nagashio, Purely in-plane ferroelectricity in monolayer SnS at room temperature, Nature Communications 11, 1 (2020).
  * Tian  _et al._ [2017] Z. Tian, C. Guo, M. Zhao, R. Li, and J. Xue, Two-Dimensional SnS: A Phosphorene Analogue with Strong In-Plane Electronic Anisotropy, ACS Nano 11, 2219 (2017).
  * Li  _et al._ [2017] M. Li, Y. Wu, T. Li, Y. Chen, H. Ding, Y. Lin, N. Pan, and X. Wang, Revealing anisotropy and thickness dependence of Raman spectra for SnS flakes, RSC Adv. 7, 48759 (2017).
  * Xia  _et al._ [2016] J. Xia, X. Z. Li, X. Huang, N. Mao, D. D. Zhu, L. Wang, H. Xu, and X. M. Meng, Physical vapor deposition synthesis of two-dimensional orthorhombic SnS flakes with strong angle/temperature-dependent Raman responses, Nanoscale 8, 2063 (2016).
  * Sarkar  _et al._ [2023] A. S. Sarkar, I. Konidakis, E. Gagaoudakis, G. M. Maragkakis, S. Psilodimitrakopoulos, D. Katerinopoulou, L. Sygellou, G. Deligeorgis, V. Binas, I. M. Oikonomou, P. Komninou, G. Kiriakidis, G. Kioseoglou, and E. Stratakis, Liquid Phase Isolation of SnS Monolayers with Enhanced Optoelectronic Properties, Advanced Science 10, 2201842 (2023).
  * Ronneberger  _et al._ [2020] I. Ronneberger, Z. Zanolli, M. Wuttig, and R. Mazzarello, Changes of Structure and Bonding with Thickness in Chalcogenide Thin Films 2001033, 10.1002/adma.202001033 (2020).
  * Raty  _et al._ [2021] J.-Y. Raty, C. Gatti, C.-F. Schön, and M. Wuttig, How to Identify Lone Pairs, Van der Waals Gaps, and Metavalent Bonding Using Charge and Pair Density Methods: From Elemental Chalcogens to Lead Chalcogenides and Phase‐Change Materials, physica status solidi (RRL) – Rapid Research Letters 15, 2000534 (2021).
  * Marini  _et al._ [2006] A. Marini, P. García-González, and A. Rubio, First-principles description of correlation effects in layered materials, Physical Review Letters 96, 136404 (2006).
  * Tritsaris  _et al._ [2013] G. A. Tritsaris, B. D. Malone, and E. Kaxiras, Optoelectronic properties of single-layer, double-layer, and bulk tin sulfide: A theoretical study, Journal of Applied Physics 113, 10.1063/1.4811455 (2013).
  * Vidal  _et al._ [2012] J. Vidal, S. Lany, M. D’Avezac, A. Zunger, A. Zakutayev, J. Francis, and J. Tate, Band-structure, optical properties, and defect physics of the photovoltaic semiconductor SnS, Applied Physics Letters 100, 10.1063/1.3675880 (2012).
  * Polo  _et al._ [2002] V. Polo, E. Kraka, and D. Cremer, Some thoughts about the stability and reliability of commonly used exchange-correlation functionals - Coverage of dynamic and nondynamic correlation effects, Theoretical Chemistry Accounts 107, 291 (2002).
  * Gonze  _et al._ [1995] X. Gonze, P. Ghosez, and R. W. Godby, Density-polarization functional theory of the response of a periodic insulating solid to an electric field, Physical Review Letters 74, 4035 (1995).
  * Gho [1997] Long-wavelength behavior of the exchange-correlation kernel in the Kohn-Sham theory of periodic systems, Physical Review B 56, 12811 (1997).
  * Perdew and Zunger [1981] J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Physical Review B 23, 5048 (1981).
  * Slater [1951] J. C. Slater, A simplification of the Hartree-Fock method, Physical Review 81, 385 (1951).
  * Perdew and Wang [1992] J. P. Perdew and Y. Wang, Accurate and simple analytic representation of the electron-gas correlation energy, Physical Review B 45, 13 244 (1992).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Physical Review Letters 77, 3865 (1996).
  * Perdew  _et al._ [2008] J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou, and K. Burke, Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces, Physical Review Letters 100, 136406 (2008), arXiv:0711.0156 .
  * Troullier and Martins [1991] N. Troullier and J. L. Martins, Efficient pseudopotentials for plane-wave calculations, Physical Review B 43, 1993 (1991).
  * Rappe  _et al._ [1990] A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Optimized pseudopotentials, Physical Review B 41, 1227 (1990).
  * Rappe  _et al._ [1991] A. M. Rappe, K. M. Rabe, E. Kaxiras, and J. D. Joannopoulos, Erratum: Optimized pseudopotentials [phys. rev. b 41, 1227 (1990)], Physical Review B 44, 13175 (1991).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Physical Review B 59, 1758 (1999).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, Journal of Physics: Condensed Matter 21, 395502 (2009), arXiv:0906.2569 .
  * Hamann [2013] D. R. Hamann, Optimized norm-conserving Vanderbilt pseudopotentials, Physical Review B 88, 085117 (2013).
  * Wan [2012] Ferroelectric Instability Under Screened Coulomb Interactions, Physical Review Letters 109, 247601 (2012).
  * Cochran [1960] W. Cochran, Crystal stability and the theory of ferroelectricity, Advances in Physics 9, 387 (1960).
  * Fei  _et al._ [2016] R. Fei, W. Kang, and L. Yang, Ferroelectricity and Phase Transitions in Monolayer Group-IV Monochalcogenides, Physical Review Letters 117, 097601 (2016).
  * Otani and Sugino [2006] M. Otani and O. Sugino, First-principles calculations of charged surfaces and interfaces: A plane-wave nonrepeated slab approach, Physical Review B 73, 115407 (2006).
  * Sohier  _et al._ [2017] T. Sohier, M. Calandra, and F. Mauri, Density functional perturbation theory for gated two-dimensional heterostructures: Theoretical developments and application to flexural phonons in graphene, Physical Review B 96, 075448 (2017), arXiv:1705.04973 .
  * Tinte  _et al._ [2003a] S. Tinte, J. Íñiguez, K. M. Rabe, and D. Vanderbilt, Quantitative analysis of the first-principles effective hamiltonian approach to ferroelectric perovskites, Phys. Rev. B 67, 064106 (2003a).
  * Tinte  _et al._ [2003b] S. Tinte, K. M. Rabe, and D. Vanderbilt, Anomalous enhancement of tetragonality in ${\mathrm{pbtio}}_{3}$ induced by negative pressure, Phys. Rev. B 68, 144105 (2003b).
  * De Silva and Corminboeuf [2014] P. De Silva and C. Corminboeuf, Simultaneous visualization of covalent and noncovalent interactions using regions of density overlap, Journal of Chemical Theory and Computation 10, 3745 (2014).
  * Damle  _et al._ [2014] A. Damle, L. Lin, and L. Ying, Compressed representation of kohn-sham orbitals via selected columns of the density matrix., Journal of chemical theory and computation 11 4, 1463 (2014).
  * Rehn  _et al._ [2015] J. Rehn, A. Sen, A. Andreanov, K. Damle, R. Moessner, and A. Scardicchio, Random Coulomb antiferromagnets: From diluted spin liquids to Euclidean random matrices, Physical Review B - Condensed Matter and Materials Physics 92, 1 (2015), arXiv:1505.01764 .
  * Fan  _et al._ [2021] Z. Fan, J. H. Garcia, A. W. Cummings, J. E. Barrios-Vargas, M. Panhans, A. Harju, F. Ortmann, and S. Roche, Linear scaling quantum transport methodologies, Physics Reports 903, 1 (2021), arXiv:1811.07387 .
  * Merkel  _et al._ [2023] K. Merkel, J. Greiner, and F. Ortmann, Understanding the electronic pi-system of 2D covalent organic frameworks with Wannier functions, Scientific Reports 13, 1 (2023).
  * Dick and Overhauser [1958] B. G. Dick and A. W. Overhauser, Theory of the dielectric constants of alkali halide crystals, Phys. Rev. 112, 90 (1958).
  * Ghosez  _et al._ [1998] P. Ghosez, J.-P. Michenaud, and X. Gonze, Dynamical atomic charges: The case of ABO3 compounds, Physical Review B 58, 6224 (1998).
  * Vanderbilt [2000] D. Vanderbilt, Berry-phase theory of proper piezoelectric response, Journal of Physics and Chemistry of Solids 61, 147 (2000), 9903137 [cond-mat] .
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, _et al._ , Advanced capabilities for materials modelling with Quantum ESPRESSO, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Umari and Pasquarello [2002] P. Umari and A. Pasquarello, Ab initio Molecular Dynamics in a Finite Homogeneous Electric Field, Physical Review Letters 89, 1 (2002).
  * Li  _et al._ [2016] L. Li, J. C. Snyder, I. M. Pelaschier, J. Huang, U.-N. Niranjan, P. Duncan, M. Rupp, K.-R. Müller, and K. Burke, Understanding machine-learned density functionals, International Journal of Quantum Chemistry 116, 819 (2016), arXiv:1404.1333 .
  * Brockherde  _et al._ [2017] F. Brockherde, L. Vogt, L. Li, M. E. Tuckerman, K. Burke, and K. R. Müller, Bypassing the Kohn-Sham equations with machine learning, Nature Communications 8, 1 (2017), arXiv:1609.02815 .
  * Louie  _et al._ [1982] S. G. Louie, S. Froyen, and M. L. Cohen, Nonlinear ionic pseudopotentials in spin-density-functional calculations, Physical Review B 26, 1738 (1982).
  * Dal Corso [2014] A. Dal Corso, Pseudopotentials periodic table: From H to Pu, Computational Materials Science 95, 337 (2014).
  * Hamann [2017] D. R. Hamann, Erratum: Optimized norm-conserving Vanderbilt pseudopotentials [phys. rev. b 88, 085117 (2013)], Physical Review B 95, 239906 (2017).
  * Fuchs and Scheffler [1999] M. Fuchs and M. Scheffler, Ab initio pseudopotentials for electronic structure calculations of poly-atomic systems using density-functional theory, Computer Physics Communications 119, 67 (1999).
  * Grimme [2006] S. Grimme, Semiempirical GGA-type density functional constructed with a long-range dispersion correction, Journal of Computational Chemistry 27, 1787 (2006).
  * Huber  _et al._ [2020] S. P. Huber, S. Zoupanos, M. Uhrin, L. Talirz, L. Kahle, R. Häuselmann, D. Gresch, T. Müller, A. V. Yakutovich, C. W. Andersen, F. F. Ramirez, C. S. Adorf, F. Gargiulo, S. Kumbhar, E. Passaro, C. Johnston, A. Merkys, A. Cepellotti, N. Mounet, N. Marzari, B. Kozinsky, and G. Pizzi, AiiDA 1.0, a scalable computational infrastructure for automated reproducible workflows and data provenance, Scientific Data 7, 1 (2020).

## Supplementary Information

###  VI.1 PDOS

(a)

(b)

(c)

Figure 10: Partial Density of states for the (a) Sn(s), (b) Sn(p) and (c) S(p)
orbitals for the different pseudopotentials used in this work (see methods
section in the main text for shorthand notation). Except for the SL, all
pseudopotentials with PBE functional are almost overlapping, and similarly for
the PBEsol and LDA approximations.

###  VI.2 Wannierization

Figure 11 shows the projectabilities of the states for SnS, as derived for PBE
SL, using the aiida software. The complementary error function (CEF) is,

$f\left(e\right)=\frac{1}{2}\operatorname{erfc}\left(\frac{\epsilon-\mu}{\sigma}\right)$
(8)

The states for which the projectability is close to one are well-represented.
Higher states (i.e. states with many radial nodes) are less easy to have a
well-localized wannier representation. The wannierization with these
parameters resulted in some localized states being outside the home unit cell
and not related to the lone pair location. Therefore, these basis states were
not used in this work.

Figure 11: Projectabilities of states as a function of the energy
(Nk$\times$Nbands points) for the PBE SL PP. The dashed line represents the
CEF. $\mu$ and $\sigma$ correspond to one choice of optimum parameter that can
be used for the CEF.