Inconsistencies in NICS/σ values of an adsorbed molecule

[jakob_schramm]

Dear developers,

I am interested in the adsorption of polycyclic conjugated hydrocarbons on metal surfaces and the influence of this interaction on their aromaticity.
I was therefore pleased to see that, starting with version 6.6.0, it is possible to calculate NICS values using periodic boundary conditions.
I subsequently began testing how this method could be applied to my research questions.
For molecules in the gas phase, the approach appears very promising, as I am able to reproduce values obtained with molecular quantum chemistry codes.
However, I quickly noted that the calculation of NICS values (or σ, because NICS = -σ) of an adsorbed molecule is highly sensitive to the chosen smearing parameters.
(Note that, in general, SIGMA needs to be small and ISMEAR cannot be >0 in order to converge the NMR calculation)
Particularly, I do not observe meaningful convergence upon decreasing SIGMA for different ISMEAR settings. Instead, the resulting values appear rather arbitrary.
Some values as an example: σ_zz(0) values of Azulene, 2 Å above a Cu(111) surface
- ISMEAR=0 & SIGMA=0.01 : -83.3 & -105.9 ppm
- ISMEAR=0 & SIGMA=0.001 : -9.5 & 50.8 ppm
- ISMEAR=-1 & SIGMA=0.01 : 38.7 & -1.8 ppm
- ISMEAR=-1 & SIGMA=0.001 : -154.8 & -73.9 ppm

Furthermore, I observed in additional test calculations that NICS values (or σ) above "diagonal" and "parallel" bridge positions are not equivalent, although they should be identical according to the threefold symmetry of the (111) surface.
Following that observation, I repeated the same Azulene adsorption calculation as before, but rotated the unit cell by 60°.
This changes the orientation of the surface and molecule with respect to the x- and y-axes, while leaving the local adsorption geometry unchanged.
I had assumed that the GIPAW approach ensures gauge independence in NMR calculations, so that such a rotation of the unit cell should not affect the results.
Therefore, the σ_zz(0) values would be expected to remain identical. But they are not:
- ISMEAR=0 & SIGMA=0.01 : -83.3 & -105.9 ppm
- ISMEAR=0 & SIGMA=0.01 & rotated cell : -57.8 & -106.3 ppm

Long story short - I have two main questions:
1. Is there any experience with calculating NICS values (or σ) for molecules on metal surfaces, in particular regarding the appropriate choice of smearing parameters? Is this meaningful from a methodological standpoint at all?
2. Why are there different NICS values (or σ) upon rotating the cell? Should these quantities not be invariant under such transformations?

Thanks a lot!
Jakob

[michael_wolloch]

Dear Jakob,

Thanks for reaching out and taking so much care to provide all the context and necessary information.

I am not an expert in NMR and NICS, but I know that NICS on metallic systems makes very little sense in VASP's implementation. I suspect that this is the reason why your results are so strongly dependent on SIGMA, and I don't think that VASP can help you with your study if metal surfaces are a prerequisite.

I will reach out to some experts at Radboud University who developed this feature to confirm my suspicion. I will get back to you as soon as I know more,
Cheers, Michael

[michael_wolloch]

Dear Jakob,

I have a somewhat more detailed answer for you after consulting with Gilles de Wijs (emphasis by me, and some slight changes to his email):

NICS is an option with shielding for NMR, so is LCHIMAG.
And LCHIMAG is only supported for non-metals, i.e., systems
with gap and all occupied states with occupancy either strictly
2 (so 2.0000 or 1.0000 for ISPIN=2) and empty states 0.00000.
So in practice with Gaussian smearing and SIGMA well below the
band gap.

Because of the presence of the Cu surface there is a metal and
one needs to consider (1) Knight shifts and (2) repair trouble
because of a vanishing denominator in the linear response.
That's not implemented. It's on the wish list.

Of course, results should be independent of system rotation.
For non-metals they are. But then be careful to set
"LNMR_SYM_RED = .TRUE." when running with symmetry,
because k-point stars (the center points of the star
are on the usual mesh) should fully map onto (other)
stars under the space group symmetry operations.
Whether this happens depends on the orientation
of the stars (that are oriented along the Cartesian
coordinate's directions) relative to the lattice vectors.
If there is trouble, the effective
IBZ is enlarged. I don't know whether this explains question 2, but the calculated shielding of a metallic surface
is anyway meaningless and ultra-sensitive to any detail.

So the simple answer: Shieldings for metallic systems
are not implemented, so the NICS in the vicinity of a
metallic surface are not meaningful.

So the message remains the same. What you are trying to do is unfortunately not possible right now with VASP 6.6.0. However, it is on the wish-list of new NMR features, and please keep a lookout for this in the changelog.

Let me know if this answers your question, and I can lock the topic.
Cheers, Michael